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ABSTRACT 

This report concerns the prediction of elastic moduli for fabric reinforced 
composites. Many analysis methods previously applied to this problem are not 
general enough for use with the more complex reinforcing geometries that have 
been suggested for suppressing impact damage and general delamination. The 
proposed analysis places no restrictions on fabric microgeometry except that 
it be determinate within some repeating rectangular pattern. No assumptions 
are made regarding fiber cross-sectional shapes or fiber paths. 

The analysis is based on a mechanical model that consists of a single 
rectangular unit cell which typically contains one repeating pattern of the 
fabric design. Every unit cell is assumed to be surrounded on six sides by 
identical unit cells, similarly oriented, with all common surfaces perfectly 
bonded. Elastic analysis of one such cell element yields all 3-D moduli of 
the composite. 

For analysis purposes the unit cell is assumed to be divisible into small 
rectangular subcells in which the reinforcing geometries are easier to 
visualize, define, and analyze. Details of micromodeling and the use of 
single, finite element hexahedra to mathematically represent each subcell are 
described in this report. The analysis is applied to a variety of woven, 
braided, and knitted fabric reinforcements in epoxy matrix. A special purpose 
computer program based on this analysis is also described. Moduli predictions 
from this program are compared to fabric reinforced composite test data. A 
program listing with sample input and output is included in the Appendices. 
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1 . I N TRQDUC TI OM 


Complex fibrous preforms (made from weaves, braids, knits, XYZ construction, etc.) 
are under consideration for use as composite reinforcing materials. As 
reinforcers, these preforms should increase impact resistance, reduce lay up time, 
and generally improve interlaminar properties. Use of these textile structures 
introduces a variety of both new and old manufacturing processes into the 
fabrication cycle and creates an intermediate material stage. This stage consists 
of an array of systematically interlaced fiber bundles (yarns or tows) surrounded 
and impregnated by a matrix material. The basic repeating structural element in 
this array has characteristic dimensions that are much larger than a single fiber 
diameter, but are on the order of the structural plate or shell wall thickness 
(Figure 1) . From an analyst's point of view this stage creates a gap between 
micro-mechanics and laminate analysis, and requires a new level of mechanics 
analysis. The term "fabric mechanics" seems appropriate for categorizing this 
level of composite mechanics. However, use of this terminology may confuse 
textile engineers who normally associate fabric mechanics with unimpregnated 
fabric behavior. The goals of fabric mechanics is to predict material properties 
at the fabric reinforced composite level from knowledge of unidirectional 
composite and bulk matrix material properties. This level of analysis is 
necessary to reduce dependence on testing and to assist in data extrapolation and 
interpretation. The influence of factors such as fiber curvature, misalignment, 
and ply nesting are not conveniently assessed by either micro-mechanics or 
laminate analysis. Fabric mechanics' role in the scheme of composite mechanics is 
shown in Figure 1. For laminates made from unidirectional material, fabric 
mechanics is not required because micro-mechanics provides all the necessary 
information for laminate analysis. 
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Although methods of fabric analysis are not firmly established, its goals are 
identical to those of all mechanics analyses; namely, to predict stresses and 
deformations within a deterministic structure as a result of prescribed external 
loads or deformations. From these predictions (plus suitable failure criteria) 
composite static properties of stiffness and strength can be estimated. The 
limited objective of this report is to develop adequate capability (and a computer 
code) to predict static stiffness of various fabric reinforced composites proposed 
for minimizing impact damage. Several methods for predicting these properties 
have been proposed, but they are built on many assumptions, and do not enjoy 
widespread use or confidence. 

The principal barrier to improved analyses is the geometric complexity of most 
fabric microstructures. The proposed analysis crosses this barrier via two stages 
of mechanical modeling. For example, consider the rib weave fabric microstructure 
shown in Figure 2. This fabric geometry can be subdivided into small block 
substructures or subcells (as shown in the same Figure) . Simpler reinforcing 
fiber geometry within each different subcell can be modeled by finite element 
hexahedra. After reassembling the subcells the boundary conditions required to 
simulate the six independent unit strain cases of 3-D elasticity may be enforced 
(as in micro-mechanics) . The resulting elastic solutions yield all necessary 
information for computation of the 6X6 stiffness or flexibility matrix of the 
composite. The approach is well-suited to computer analysis. A special propose 
code has been written that predicts the stiffness of fabric reinforced composites 
on the foregoing basis. A laminate analysis must be performed if several 
different fabrics are combined into a thick, laminated structure. This approach 
has the added advantage of lying within the range of experience of most structural 
analysts . 
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One of the long range objectives of this work is to develop the capability of 
estimating all of the 3-D static material property inputs required of transient 
dynamic impact codes that analyze the response of composite target plates to 
spherical impactors (Ref. 1). This capability will permit rational selection of 
woven, braided, and stitched fabric parameters to best meet impact damage 
requirements with the least penalty to the basic static design properties. 
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2. REINFORCING GEOMETRY 

In the textile sense, the term "fabric" is applicable to any network of bonded or 
unbonded fibers which results in an essentially planar structure that is flexible 
yet cohesive beyond the point of being self-supporting. The fabrics can be 
categorized as having either a random distribution of fiber (as in a felt or 
bonded flock) or an orderly distribution (as in weaves, braids or knits) . Hybrids 
of both are possible and some elements of order and disorder are always present in 
both. Yarns or tows can also be continuous or discontinuous or a hybrid of both 
(as in cut pile carpets) . In this report, only fabrics made from continuous tows 
arranged in identifiable and repetitious patterns of entanglement will be 
considered. It is also assumed that this repetitious pattern has finite 
dimensions and can be contained exclusively within the volume of a rectangular 
prism. This prism is termed a "unit cell" of the fabric. Reproduction of this 
unit cell, rigid transformation of it, and attachment to compatible surfaces of 
other unit cells, allows any large planar area to be covered by this fabric, with 
areal contours closely approximated by outer boundaries of units cells. Also unit 
cells can by layered to fill any closed volume, even where minimum dimensions of 
that volume are several orders of magnitude larger than unit cell dimensions. It 
is assumed in this process that all common boundaries between unit cells are 
perfectly bonded to each other. 

This 3-D unit cell often visually relates to the 2-D fabric design plan that, 
along with the harness draft and the chain draft, is typically used to describe a 
weave pattern and to set up a shuttle loom to make that weave. The design plan is 
a coded representation of a weave, whereas the unit cell is a true scale model of 
the weave, albeit a somewhat idealized one. A design plan is usually drawn on 8 X 
8 design paper. Columns represent warp yarns; rows represent fill yarns. Each 
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small square represents an Intersection of one warp and one fill yarn. A shaded 
square indicates that the fill yarn crosses over the warp yarn at that 
intersection. 

One unit cell representation of a plain weave is shown in Figure 3A. The 
corresponding fabric design plan is shown in Figure 3B. But, neither the unit 
cell nor the design plan are unique. Another very different unit cell of the same 
weave is shown in Figure 3C. 

Since unit cells are not unique, the question arises as to which unit cell should 
form the basis for fabric analysis. The answer is that unit cell which is 
simplest to analyze. Often, the simplest cell is presumed to be the smallest 
cell. However, the smallest unit cell is frequently difficult to analyze, 
particularly when it has complex boundaries, for example the twill weave (Figure 
4) . Boundaries of the smallest twill weave unit cell do not form a rectangular 
prism. Hence, mechanical analysis of this weave can be difficult. The larger 
unit cell of Figure 4C does not pose this problem. 

Several rectangular unit cells for common weaves, knits and braids are shown in 
Figure 5. For some complex fabrics, it is difficult to identify a rectangular 
unit cell. For woven fabrics, one repeat of the design plan often resembles an 
out-of-scale projection of a unit cell surface. But with knits or braids there is 
no corresponding 2-D design plan. However, rectangular unit cells generally can 
be established, although some liberties with the tow geometry may be required. 
For example, some tow geometries like satin weaves (with large, complex, 
rectangular, unit cells) contain smaller, simpler, trapezoidal, unit cells. 
Slightly distorting the boundary of the smaller unit cell can make it rectangular 
but this process entails some violation of the true microgeometry. 
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3. HISTORICAL DEVELQPMEN!£ 

Woven fiberglass cloth was the first reinforcing material with enough strength and 
stiffness for general use on aircraft structures. However, relatively few fabric 
weave styles saw extensive application. These fabrics were experimentally 
characterized using epoxy and polyester resins and were treated like new 
orthotropic materials with little consideration of fabric microstructural 
analysis. Although it was known at that time that fabric microgeometry had a 
significant effect on moduli and was one of the determining factors in stress 
failure of this material, a convincing mechanical analysis of this problem was 
beyond computational capabilities. The advent of filament winding and 
unidirectional prepregging, along with increased capability in computing, soon 
made such an analysis possible at the unidirectional level. It was assumed that a 
wide range of fiber volume fractions and fiber cross-sectional shapes would 
necessitate some type of micro-mechanics design/analysis cycle. Instead 
unidirectional materials became even more standardized than fabrics. Test data 
replaced micro-mechanics. 

The search for improved impact resistance has brought about a renewal of interest 
in more complex reinforcing geometries. Meanwhile, computing capability has kept 
pace with materials development. This increased computing capability now allows 
design and analysis of fabrics at a level comparable to design and analysis of 
unidirectional materials. The number of geometric variables present and the 
requirements of impact resistance suggest fabric reinforcements will play a role 
in developing impact resistant composites. 

First attempts to extend composite mechanics beyond the realm of planar analysis 
arose in conjunction with carbon-carbon development for reentry vehicles (Ref 2). 
These attempts ranged from 3-D versions of the rules of mixture and netting 
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analysis to finite element solutions. Fiber curvature was seldom considered. 

Two-dimensional analysis of wavy fiber reinforced composites was considered by 
some to be an initial step toward the analysis of true fabric reinforced 
composites (Ref. 3) . In this approach a single wavy fiber is treated as a curved 
beam on an elastic foundation. 

Another step toward geometric reality (Ref. 4) involved modeling of orthogonal 
fiber arrays. However, the problem was reduced to two dimensions and solved by 
finite elements. 

Among the first to treat fabric mechanics as a true 3-D problem were Ishikawa, 
Chou, et al., (Refs. 5 through 11). Based on assumed fiber paths and cross- 
sectional shapes these authors developed a variety of models based largely on 
"strength of materials" type idealizations and extended some of this work into the 
nonlinear property realm. 

Materials Sciences Corporation developed a simple fabric model based on earlier 
work on random fiber reinforcements (Ref. 12) . 

Another more refined approach was the attempt to fit fabric mechanics directly 
into the context of laminated plate theory (Ref. 13) . This analysis also 
occasionally relies on "strength of materials" assumptions (Ref. 14). 

The primary legacy of the previous work lies not in analytical details, but in the 
modeling of complex fabric microgeometry with unit cells and subcells. Not only 
does this reduce the geometry problem to smaller, "more digestible" pieces; it 
also meshes with the finite element scheme of reducing a complex structure to an 
assemblage of simpler ones which can be characterized, reassembled, and 
manipulated to predict the response of the unit cell and the composite. 
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4. FABRIC ANALYSIS 

For convenienee in establishing thickness direction properties, each fabric 
reinforced composite ply hereafter will be considered to be imbedded in a thick 
laminate of identical plies layered such that each unit cell in each ply 
represents one level in a vertical stacking of unit cells extending through the 
full thickness of the laminate (like stories in a multi-story building) . For 
example, treat each unit cell as a single brick in a stack bond wall (Figure 6) . 

A wide variety of other stacking arrangements are equally valid and can be 
observed in laminates. Cursory investigations show these stacking variations to 
be of secondary importance to the principal ply stiffness. A more thorough 
investigation would be appropriate to strength analysis but of lower priority to 
this study. 

— Determination E lastic Constants 

Now consider a large cube of fabric reinforced composite material in which 
dimensions of the cube are orders of magnitude greater than dimensions of the unit 
cell. For most engineering purposes such a composite material can be considered 
homogeneous and anisotropic with a linear stress-strain law of the form; 
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where: 
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ii 



a 

ii 


X 

ij 


S 

ij 


= average normal strain, 

= average shear strain, 

= average normal stress (psi) 
average shear stress (psi) 

2 

= flexibility coefficient (in /lb) . 


If the cube of material is subjected to a unit average strain in the x- 
direction, with all other average s\rains held to zero, then CT^x*’ 

Xy 2 / "^xy represent average stresses corresponding to this strain state. In a 


hypothetical test these stresses could be considered to be measured quantities. 


The stress-strain law given above would then consist of six independent equations 
in the 36 unknown coefficients. Repeating this test with ^yy ^ ^ 

other average strains held at zero would yield six more equations in the . 


unknowns. All six independent unit strain states would yield 36 equations in the 

36 unknown S... These terms are the flexibility coefficients of the composite. 
13 

Inverting the matrix of flexibility coefficients yields the stiffness matrix [D] ; 
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where: D. . = stiffness coefficient of composite (psi) . 

13 

E , V. , G. . etc. for the composite can be obtained by considering the material 
ii X3 X] 


response to single nonzero components of each of the six average stresses 


where : 


E.. = modulus of elasticity of material (psi), 

11 

V,, = Poisson's Ratio, 

11 


G. . 
11 


= modulus of rigidity of material (psi) . 


4.2 Unit Cell BoundarV-Conditlons 


Now return to the initial problem of obtaining the average composite stress 
components CTxx, ^yy, ^zz, ^yz, "^xz, '^xy resulting from the average strain state 


E =1, e = e = Y = Y « Y = 0. 

XX yy zz * yz * xz * xy 


(4.3) 


This problem can be resolved by considering a single unit cell of the material, 
sufficiently far removed from the surfaces to be free from edge, corner, or 
surface effects. Since all such unit cells in the material are indistinguishable 
from one another the response of any two unit cells to any uniform average stress 
or strain state must be similar, except for rigid body motion. 


Relative displacement of the eight corners of the unit cell are shown in Figure 7. 
Any other set of relative displacements would be incompatible, inhomogeneous (in 
the average strain sense), or in violation of the unit strain case that it 
represents. 


Now consider boundary conditions on surfaces of the unit cell. First consider 
interlamina upper and lower surfaces of the cell. For every point on the upper 
surface there is a corresponding point directly below it on the bottom surface 
(Figure 8) . Call these two points "image points" because the top image point on 
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one unit cell is coincident with the bottom image point on the unit cell directly 
above it. Similarly the bottom image point on any cell coincides with the top 
image point on the unit cell below it. Since there is no reason to expect any two 
cells to deform differently under uniform average strain, the displacements of 
image points can differ only by the amount of displacement associated with the 
homogeneous unit strain case (plus some arbitrary rigid motion). Ignoring rigid 
motion, the displacement boundary conditions (corresponding to unit strain case 
4,3) for image points i,j (Figure 8) have the form: 


U. - U., V. = V., W. = W. 

131313 


(4.4) 


where U,V,W designate displacements in the positive x,y,z coordinate directions, 
respectively. This condition provides three equations in the six displacement 
components of the two image points. Three additional equations are required for a 
complete statement of a 3-D elastic problem. 


Stress continuity considerations provide these equations. The stress vector at any 
point on the lower surface of one unit cell is the reaction to the stress vector 
at the same point on the upper surface of the unit cell lying directly below it. 
Also, surface stresses on one unit cell are assumed to be identical to stresses on 
any other cell. Therefore, stress vectors at image points of the same cell must 
be equal in magnitude but oppositely directed. These conditions 
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where CT,X represent surface stress components, provide three more equations at 
each pair of image points. These conditions (4.4 and 4,5) place the problem in 
the category of a mixed boundary value problem in elasticity; i.e. both stress and 
displacement conditions apply. A similar argument may be applied to image points 
on any other pair of opposite sides of the unit cell. On sides perpendicular to 
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the x-axis the U displacement condition has the more general form Uj + 

where A is a constant that accounts for the stretching in the x~direction. 
Similar conditions can be specified for the other two normal strain cases. 

Boundary conditions for the three unit shear strain cases differ from normal 
strain cases only in that the constant terms that appear in some of the 
displacement boundary conditions are associated with displacement components that 
are parallel rather than normal to sides of the unit cell where they apply. 

4.3 Symm etry Boundary Conditions 

Structural symmetry often can be used to reduce the portion of the unit cell that 
needs to be analyzed. However, boundary conditions on the plane of symmetry, and 
the cell face parallel to the plane of symmetry, will differ from those conditions 
previously stated. Consider only the special case where two parallel faces of the 
reduced unit cell are both planes of symmetry. 

First consider the three unit normal strain problems. All eight corners of the 
reduced (by symmetry) unit cell now lie on planes of structural (and loading) 
symmetry, but their pure displacement boundary conditions are not changed from the 
previous discussion; i.e. their displacements are what establish the particular 
unit strain case under consideration. Points lying on one plane of symmetry have 
the traditional symmetry conditions of zero (or constant) normal displacement and 
zero shear stress. Thus, stresses and displacements at image points that lie on 
parallel planes of symmetry are no longer related through their boundary condition 
equations. 
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Unit shear strain load cases divide into two categories: symmetric loading 
categories and antisymmetric ones. Symmetric loading cases do not differ from 
unit normal strain cases. Antisymmetric loading categories have boundary 

conditions that are the converse of the symmetric ones; i.e., displacements 
parallel to a plane of structural symmetry are zero (or constant) , and stresses 
normal to that plane are zero. 
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5 ■ SUBELEMENT MECHANICS 

The previous section (4.) of this report describes the elastic structure (the unit 
cell)/ the various load (or displacement) cases applied to that structure, and the 
boundary conditions prevailing on surfaces of the structure. This section 
considers a numerical method for estimating internal (and surface) displacements 
of the unit cell for each set of loads and boundary conditions. This analysis is 
a direct application of the finite element method in which the unit cell 
structure is idealized into a number of smaller polyhedra (subcells) joined 
together at common vertices. The presumption is that stress, strain, and 
displacement fields are simpler to approximate within the subcell, than in the 
larger structure. Each subcell is then analyzed for the most general set of loads 
(or displacements) at each of the vertices and reassembled into the unit cell 
structure. 

If the unit cell is in the shape of a rectangular prism it is convenient to 

consider only subdivisions of the unit cell into smaller rectangular prisms 
(subcells). The advantage of this restriction is that it simplifies the geometry 
and the analysis, ' and eases coding for digital computation. This restriction is 
one of convenience and is not particularly restrictive because fabric 

microgeometries are often rectangular to some degree. For example, consider the 
plain weave unit cell of Figure 9.* This unit cell of Figure 9 can be subdivided 

into 16 subcells as shown. These subcells are all rectangular prisms, like the 

unit cell. 

* A subsequent section contains photomicrographs of actual fabric reinforced composites, induding a plain weave. The 
similarities between cross sections of the more realistic unit cells which model the composite microstructure and the 
(^)liotomicrographs are apparent. 
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If the warp and fill tows are similar and their paths are the same, then the 16 
subcells are all simple transformations of each other; i.e., any one of their 
reinforcing microgeometries can be reproduced by a combination of rigid motions 
and reflections about the face planes of any one of the subcells. The 
transformation makes the chore of describing microgeometry much simpler since a 
detailed description of one subcell (mastercell) and a set of instructions for the 
series of coordinate transformations will suffice to reconstruct the unit cell. 
As another example consider the unit cell of the 2/2 rib weave construction with 
similar warp and fill yarn geometries (Figure 2). This unit cell can be 
constructed from a sequence of transformations of two master subcells designated I 
and II in Figure 2. 

5.1 Subcell Analysis. 

Now consider the elastic analysis of a single rectangular subcell. For 
incorporation into a finite element analysis at the simplest level, it is 
necessary first to obtain a stiffness matrix relating the three components of 
force, acting at each corner of the subcell, to the three components of 
displacement at each corner. The method of obtaining this matrix is the central 
problem in fabric reinforced composite stiffness analysis. Many possible 
approaches exist. Simpler approaches often overlook fiber reinforcing detail. 
Refined analyses lead to greater computational efforts than the problem warrants. 
Hence, a compromise is needed. One compromise is suggested by the general energy 
formulation for the stiffness matrix of rectangular finite element hexahedra (Ref. 
1 5 ) ; i . e . , 
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IP 


[k] = J i i [B] tD) [B] d(vol) 

VO. 


(5.1) 


where: [D] = 3-D material stiffness matrix, 

[B] = strain/displacement matrix, 
vol = volume 

[D] contains only local material property distribution functions and [B] contains 
only derivatives of displacement mode shapes. Superscript T designates matrix 
transposition. [D] is usually obtained by inverting the flexibility matrix [S] . 
(Si is usually obtained by transforming the stresses and strains from material 
coordinates of orthotropy into the coordinate of integration. In the natural 
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(5.2) 


V(x,y,z), and W(x,y,z). A variety of displacements have been associated with 
hexahedra. The first hexahedra were formed from various combinations of 
tetrahedra whose displacements were assumed to be linear. Higher order tetrahedra 


^Assuming the constitent materials to be orthotropic rather than generally anisotropic does not restrict the response of the 
unit cell to be less than generally anisotropic 
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were also used. Later, isoparametric families of hexahedra were used with 8, 20, 
32 nodes. All nodes were located on element edges. Elements with internal nodes 
were also investigated. Also used were super-elements formed from smaller 
hexahedra. Displacement derivatives were used as nodal degrees of freedom also. 
The intent of this study was to begin with simple elements, increasing their 
complexity only as required to obtain acceptable accuracy for moduli predictions. 
The first displacements chosen were those associated with eight-node isoparametric 
hexahedra (Ref. 15); i.e.; 

8abc(U,V,W}=(a+2x) (b-2y) (c-2z){U^,V^,W^) 

+(a+2x) (b+2y) (c-2z) {U 2 ,V 2 ,W 2 > 

+(a-2x) (b+2y) (c-2z) (U^,V 3 ,W^} 

+(a-2x) (b-2y) (c-2z) 

+(a+2x) (b-2y) (c+2z) 

+(a+2x) (b+2y) (c+2z) (U^,V^,W,) 

boo 

-Ma-2x) (b+2y) (c+2z) 

+(a-2x) (b-2y) (c+2z) {Ug,Vg,Wg} 

where: 

a == subelement side length in x-direction, 

b *= subelement side length in y-direction, 

c = subelement side length in z-direction. 

Later, additional internal degrees of freedom associated with incompatible mode 
corrections (Refs. 16, 17) were incorporated into the current analysis. 

All elements of the [B] and [D] matrices are functions of spatial coordinates 
within the subcell. If more than one material is present in the subcell, [D] is 
discontinuous in the region of interest and derivatives of true displacements will 
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also be discontinuous. Nevertheless, the integration for obtaining the stiffness 
matrix can still be performed. Mathematically, functions of class (without 

continuous first derivatives), can be approximated by a series of functions of 
class (with continuous first derivatives) . 

Formation of the stiffness matrix from [B] and [D] requires numerical integration 

T 

of the matrix product [B] [D] [B] . Integration of a matrix product is the 

integration of each element of the product matrix. Various mensuration formulae 
exist which are applicable to the formation of this integral; e.g., Gaussian or 
Newton-Cotes quadrature. In this report only the simplest step-function 
approximation to the integrand is used for computing stiffness matrices. The 
integral is formed independently for each different material appearing in the 

subcell, and then is summed over all materials in the subcell. 

For integration purposes, a 3-D grid is superimposed on the subcell volume as 

follows. First, each set of parallel subcell edges is divided into N unequal 

segments. N can be different for each coordinate direction. Planes are 

constructed through these points of subdivision, normal to the edges of the 

subcell. Intersections of these three sets of planes, within (or on) boundaries 

of the subcell, are the integration points (Figure 10) . Each element of the 

T 

[B] [Dl [B] matrix is evaluated at each integration point. 

For each material a local average value of a fiber reinforcing direction is input 

at each integration point. Each average value is estimated by considering all of 

the material closest to that integration point. The fiber direction is 
established by specifying the two spherical angles that the fiber 

principal axis makes within a local coordinate system parallel to the axes of 
integration (See Figure 11) . Elastic properties of the material are then 
transformed into the subcell reference system of integration. The value of each 
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element of the matrix product [B] [D] [B] , for that material, is then computed at 
the integration point and multiplied by the product of the volume fraction of the 
material at that point times the total rectangular volume of neighboring space 
associated with that integration point (Figure 12A) . The resulting matrix 

quantity is then totaled over each integration point and material. 

Since the integration grid spacing is variable, the volume of any material 

associated with any integration point can be obtained by summing the material 
volumes of the eight octants of volume adjacent to the integration point (See 
Figure 12B) . 

For simplicity it is often convenient to consider one volume of impregnated tow 
(for example a warp tow) as one material, while another volume of impregnated tow 
(for example a fill tow) is considered to be a different material, even though 
they may have the same elastic properties. This avoids the assigning of more than 
one fiber direction to a material at an integration point. Many details of 

reinforcing microgeometry are lost in the integration process. Examples include 

the fiber location with respect to an integration point and local fiber curvature. 

5.2 Subcell Transformations 

As mentioned previously, the ability to transform the stiffness matrix of a 

subcell into the global coordinate system of the unit cell is essential. For 

rectangular hexahedral subcells, only a limited number of transformations are 

needed, namely, rigid rotations about the coordinate axes and reflections about 

those axes. First consider a reflection about the x-coordinate axis of Figure 13. 

This is the coordinate transformation x = -x, y “ y, z = "z. Let the node points 
of the subcell be ordered in the sequence shown in Figure 13 and forces (X^,Y^,Z^) 

and displacements associated with those nodes ordered as follows: 
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{Uj,Vi,Wi,U2,V2,W2, ••• ^«8'''8'”8^* 

The rearrangement of the microgeometry caused by reversal of the x-axis is 
reflected in the rearrangement of rows and columns of the [k] matrix as follows: 


{ X^, -Xg, Yg, Z^, ~X^,Y^,Z^, Yg, Zg, 


{-U5,V3,W5,-Ug,Vg,W6,-U7,V^,W^,-U8,V8.W8, 


-U^,V^,W^,-U2^V2,W2,-U3,V3,W3,-U^,V^,W^} 


Similar exchanges of rows, columns, and signs can characterize a reversal of fiber 
geometry about the y or z-axes. 

Now consider changes in fiber reinforcing geometry brought about by rotation of 
the subelement about a coordinate axis. Look at rotation about the x-axis of 90^ 
as shown in Figure 14. The following rearrangeirient of rows and columns of the [k] 
matrix achieves this reinforcing geometry modification: 


{X^,Z^,-Y^,X^,Z^,-Y^,X2,Z2,-Y2/X3,Z3, Y3, 
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A similar rotation of 180 about the x-axis leads to the following rearrangement 
of forces and displacements: 


(U3, -V3, -W3,U^, -V^, -W^, U3, -V3, -W3, U3, -V3, -w^, 

«7 ' -'^7 ' -”7 ' “8 ' -'^8 ' -"8 ' “5' -'^S' -"5 ' "6 ' -'^6' '"e ' 


A rotation of 270 about the x-axis is produced by the following row/column 
substitution: 


iU3,-W2,V3,U3,-W3,V3,0^,-W^,V^,Ui,-Wi,V3, 

“6'-”6''^6'"7'-”7'^'''8''”8'''8'“5'-"5'''5> 


In an analogous manner, stiffness matrices corresponding to subelement rotations 
or reflections about the y and z-axes may be obtained. Any sequence of rotations 
and/or reflections may follow one another in order. 
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6. ASS£MB.LI .. . AN.D. .^QLUXXQM 

The process of assembling various subcell stiffness matrices into a unit cell 
stiffness matrix is similar to any other finite element assembly. First, nodal 
forces and displacements are arranged in some convenient order. Each element of 
each subcell stiffness matrix is then placed in its appropriate location within 
the unit cell stiffness matrix. All terms in any location of the unit cell matrix 
are summed to obtain the unconstrained stiffness matrix of the assembled unit 
cell. Now consider the manner in which the nodal displacements are obtained from 
the unit cell stiffness matrix, the applied displacements (at corners of the 
cell), and the surface boundary conditions. 

6.1 Discrete Unit Cell Boundary Conditions 

Unit cell corner displacements are established by the particular applied unit 
strain and by the location of the zero displacement reference axis. In this 
report the unit cell is always located wholly within the positive octant of a 
right hand x,y, z coordinate system with one corner of the unit cell at the origin 
(Figure 7), the zero displacement reference point. If the unit cell has 
dimensions A,B,C in the x,y, z directions respectively, then the eight corner nodes 
have the coordinate locations shown in Figure 7. Corner node displacements for 
the six average unit strain cases of interest are given in Table 1. 

Next consider boundary conditions that apply at a finite element node point 
situated on either the upper or lower surface of the unit cell, a common point on 
two adjacent plies (Figure 8) . The nodal boundary conditions are equivalent to 
the continuum elastic boundary conditions discussed previously (Section 4.), 
Nodal forces at two image points (i,j) must be equal but oppositely directed 
(Figure 8) , i . e . , 
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X. - -X., - -Yj, - -Z^ 


( 6 . 1 ) 


while the corresponding displacement boundary conditions at the same two nodes 


are: 


Unit 

Strain 

Case 

(1) 

E 

XX 

s= 

1*; 

u. 

1 

« U.; 
J 

V. 

X 

“ V 

’'i 

= w. 
3 

\ 

Unit 

Strain 

Case 

(2) 

e 

yy 

xa 

1*; 

u. 

X 

» U.; 
3 

V. 

X 

- ''y 

w. 

X 

= W. 
3 


Unit 

Strain 

Case 

(3) 

e 

zz 


1*; 


-Uy 

V. 

X 


w. 

X 

-C = 

w. 

3 

Unit 

Strain 

Case 

(4) 

T 

* yz 


1*; 

"i 


V. 

X 

-C * 


W. = 

X 

w. 

3 

Unit 

Strain 

Case 

(5) 

y 

* xz 

= 

1*; 

u. 

.X 

-C “ 

"j'- 

V. = 

X 

V . ; 
3 

W. = 

X 

W. 

3 

Unit 

Strain 

Case 

(6) 

7 

* xy 

sss 

1*; 

u. 

X 

= U,; 
3 

V, 

X 

*= V . ; 
3 

w. 

X 

= W. 
3 

j 




( 6 . 2 ) 


If there is more than one subcell per ply thickness then there may be node points 
on all four vertical faces of the unit cell that are not edge or corner nodes. If 
these vertical faces are not planes of structural symmetry then the mixed boundary 
conditions at these node points will be similar to those just stated for upper and 
lower surface image points. Only the magnitude and location of the constants in 
the displacement constraints will differ. 

6.2 Discrete Symmetry Boundary Conditions 

Where structural symmetry conditions exist on a vertical face of the unit cell, 
nodal force conditions become, for all nodes (except corner nodes) on planes of 
symmetry perpendicular to the 

*AII other average strain components equal zero. 
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X-axis 



(6.3) 


(6.4) 


Displacement constraints, on planes of symmetry perpendicular to the x-axis, are: 

U « 0 (or A) for Unit Strain Case (1) 

U = 0 for Unit Strain Case (2), (3), (4) 

V = w * 0 for Unit Strain Case (5), (6) 

For planes of symmetry perpendicular to the y-axis: 

V = 0 (or B) for Unit Strain Case (2) 

V = 0 for Unit Strain Case (1)^ (3), (5) 

U = W “ 0 for Unit Strain Case (4) , (6) 

Similar sets of boundary conditions can be specified for all other unit average 
strain cases with the chief difference being size and location of the constant 
terms in the displacement constraints. The general form of the system of 
equations for nodal forces and displacements can now be written. 



The assembled unit cell stiffness matrix, [K] , relates all nodal forces, {F}, and 
displacements {8}, as follows: 

(F) ^ [K] (8). (6.7) 

The homogeneous force boundary conditions at surface and edge nodes can be written 
symbolically in matrix form as: 

(0) = [H] {F} (6.8) 

where (O) is the null vector and [H] the coefficient matrix. 

The homogeneous and inhomogeneous displacement boundary conditions at surface, 
edge and corner nodes can similarly be written in matrix form as: 
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<0,A,B,C} - [J] {6> 


(6.9) 


where {0,A,B,C} is a vector containing only the values 0,A,B or C and [J] is the 
coefficient matrix. Implicit is the absence of force resultants at all internal 
nodes . 

The foregoing system of equations (6. 7, 6. 8, 6. 9) has dimension 6n by 6n where n is 
the number of node points in (or on) a unit cell . The size of the system can be 
substantially reduced before attempting a solution. 

All nodal forces must be eliminated by first setting to zero those forces that are 
zero in the stiffness matrix. Remaining nonzero force equations involving the 
Stiffness matrix are then used to eliminate all nonzero force terms from the 
homogeneous force boundary conditions. The combination of zero-force stiffness 
equations and the force boundary conditions (in terms of displacements) comprise a 
set of homogeneous equations in nodal displacement variables only. Each 
displacement boundary condition is so simple that it can be used to eliminate one 
surface displacement variable in the foregoing set of homogeneous equations. Over 
half of the displacement variables on the faces of the unit cell are eliminated 
in this manner. The remaining equations can now be solved for all unconstrained 
displacements. The foregoing procedure is not standard finite element procedure. 
Most general purpose finite element programs do not treat boundary conditions of 
this complexity. 

Displacement boundary conditions and corner conditions are then used to obtain 
constrained displacement values. All surface and corner forces can be obtained by 
substituting the complete set of nodal displacements into the original stiffness 
equations. Nodal force components on any side of the unit cell can then be summed 
and divided by the area of the side. This dividend is the average stress on that 
side. These stress values are the coefficients of the 3-D composite stress/strain 
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law, [D] . 


Elastic solutions for the six unit strain cases provide all 36 
coefficients of the stress/strain law which can be inverted to provide the 
composite flexibility matrix. Recourse to definitions Of engineering constants 
provides formulae for these quantities. 

One remaining simplification eliminates the need for reformulating and resolving 
each elastic displacement problem for each of the six unit strain cases. Each 
unit strain case is considered to be a combination of a homogeneous (uniform) 
strain field plus an inhomogeneous (nonuniform) one. Nodal displacements for the 
six uniform strain fields are obvious from knowledge of the response of 
homogeneous material to uniform surface stress states. The six sets of equations 
for the nonuniform displacements turn out to differ only in the inhomogeneous 
portions of their algebraic equations. Thus only one large matrix inversion is 
required (as opposed to six) . If NX is the number of subcells along the x- 
parallel side of a unit cell and NY is the number along a y-parallel side then the 
size of the inversion is 3 [ (NX) (NY) -1 ] square. 

The use of planes of structural symmetry to reduce the amount of unit cell 
structure that must be analyzed complicates matters somewhat. Strain eases 
representing symmetric loads (with respect to planes of symmetry) do not reduce to 
the same set (or number) of unconstrained displacement equations as the shear 
strain cases which represent asymmetric loadings. Thus, the computational effort 
is approximately doubled. Despite this increased effort, it is still advantageous 
to use structural symmetry when it exists. 
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7 ■ COMPUT ER PROGRA M DESCRIPTION 

A computer code, based on the foregoing analysis, has been written in Fortran 77 
language for use on the CDC 6600 computer at NASA Langley Research Center's 
central computing facility. A listing for the program is given in Appendix I. 
The core of the program consists of a set of three nested do loops. The outer do 
loop computes the stiffness matrix for each elastically different subcell in the 
finite element model of the unit cell, transforms it into the global unit cell 
coordinates, and inserts each element of the subcell stiffness matrix into its 
proper location in the larger unit cell stiffness matrix. 

The intermediate do loop ranges over each material within the subcell. In other 
words, the subcell stiffness matrix is assembled, one material at a time, and the 
contribution that each material makes to the stiffness matrix is computed and 
added into the subcell stiffness matrix, which maintains a running total of each 
material contribution, in the same manner that the unit cell stiffness matrix 
consists, at any point in time, of a running total of the different subcell 
contributions. ’ 

The third and inner most do loop ranges over each integration point in the 
integration grid. The contribution of one material to one subcell stiffness 
matrix represents the sum of its contributions at each integration point. As the 
fiber directions and constituent material volumes are read in, for any integration 
point, the contribution to the appropriate subcell stiffness matrix is calculated 
and added into the proper matrix locations. As the last constituent material 
associated with the last integration point for the last subcell is read into the 
program the last contribution to the unit cell stiffness matrix is put into place 
and the unconstrained stiffness matrix of the unit cell is completed. There is no 
further input data required. The six unit strain problems are solved in sequence 
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and the mean stresses corresponding to each unit strain tabulated. The elastic 
constants for the composite are then computed and displayed. 

Prior to entering the nested do loops the necessary information to characterize 
the constituent materials, set the number of master subcells, size the unit cell, 
divide the unit cell into subcell compartments, and fill each compartment with a 
properly transformed subcell, must be specified. The subcells are sized in the 
outermost portion of the nested do loops. The input data details are discussed in 
Appendix II. Appendix III contains an interactive sample problem input and output 
for a simple woven composite architecture. Figure 15 contains a flow chart of the 
program. 

The math subroutines that invert matrices and solve simultaneous equations (MATOPS 
and GELIM) are not included in the listing. Similar subroutines are included in 
any math library package. 


28 



2 ^ ai^PLE applications 

In this section fabric analysis is applied to a set of problems that are not 
representative of any actual fabric reinforced composite, but are simple enough to 
illustrate some basic features of the analysis. The simplest possible application 
is the prediction of the stiffness of a bulk orthotropic material, for example, a 
graphite/epoxy with the following elastic properties: 
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( 8 . 1 ) 
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where planes perpendicular to the fiber axis are planes of transverse isotropy. 
The fiber direction is denoted by 1, the direction perpendicular to the fiber (in 
the plane of the laminate) by 2, and the ply thickness direction by 3. 


The unit cell of this reinforcing geometry can be any size of rectangular prism 
with a minimum side length several orders of magnitude larger than the average 
fiber diameter. This unit cell can be assembled from a single subcell which 
contains only one material, the graphite/epoxy composite with fibers paralleling 
the x-axis. The subcell is the unit cell. The simplest 3-D integration grid and 
integration scheme is adequate in this case. Let the eight corners of the subcell 
be the integration points. The subcell stiffness matrix generated by integration 
represents a homogeneous material uniformly distributed throughout the subcell. 
The unit cell stiffness matrix is identical to the subcell stiffness matrix. 
Subsequent calculation of the composite stiffness coefficients and E,V,G values 
yields the same numbers as the material property input. Analyzing the same 
material with the fibers at an angle to (rather than parallel to) the global x- 
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axis gives the composite moduli in a different reference system from the input 
data . 

8.1 Fibe r Path Variation 

The effect of ±45° zigzag ( MA' or rickrack fiber paths on the moduli can also 
be evaluated with this analysis. For example, consider the 2-D rickrack 
reinforcing pattern of Figure 16A, The unit cell shown contains one repeat of 
this geometry. The unit cell can be subdivided into two subcells (I and II) with 
constant fiber direction in each (Figure 16B) . One subcell can be generated from 
the other via a 180° rotation about the z-axis. Thus only one master subcell and 
a pair of rigid transformations are required to orient fibers properly in the unit 
cell. Using the previous unidirectional composite properties (Equation 8.1), the 
analysis of the unit cell gives the following results; 
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where the coefficient of mutual influence > designates shear strain in the 

, k 

ij plane as a result of a unit normal strain in the k direction. The comparable 
results for unidirectional material oriented at an angle of ±45° to the x-axis (in 
the xy-plane) without the rickrack pattern are; 
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The Young's moduli and Poisson's ratios of the rickrack pattern are the same as 
the skewed unidirectional material. However, the G„ shear modulus and all of the 

non~zero coefficients of mutual influence differ. 

Neither of these examples shows how multiple materials can be mixed in the same 
subelement. The previous example problem could have been approached in this way; 
i.e., the unit cell could have been modeled by a single subcell containing both 
fiber angles in the rickrack pattern (see Figure 16C) . If half of the subcell 
contained only + 45^ material and the other half -•45*^ material then (using the 
same eight corner integration points) the volume of material associated with the 
four integration points on the x * 0 plane would consist only of +45^ material 
while integration points on the x =» A plane would have only -45° material 
associated with them. The complete subcell stiffness matrix includes all eight 
integration points and both + 45° and -45° material. Resulting moduli predictions 
are identical to the previous model based on two homogeneous subcells per unit 
cell . 

Three-dimensional rickrack reinforcements, shown in Figure 17A, can be analyzed 
almost as easily as 2-D ones. A unit cell of this fiber reinforcing pattern is 
shown in Figure 17B. Solid lines indicate fiber direction in an entire subcell. 
All four subcells (Figure 17C) are rigid transformations of each other. Variation 
in Young's modulus of this material with the z-axis fiber orientation angle (for 
the same graphite/epoxy) is shown in Figure 18. This figure also contains a plot 
of Young's moduli for the 2-D rickrack pattern of Figure 16. Internal constraints 
from adjacent subcells increase the 3-D moduli very slightly over the 2-D case. 

8.2 Voi d Analysis 

Consider another example problem in which the material is homogeneous and 
isotropic except for a periodic array of continuous, parallel, square holes. The 
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major axis of these holes is parallel to the z-axis, as shown in Figure 19A. One 
unit cell of this microgeometry consists of a cube of material with one hole in 
its center (Figure 19B). Since there is no variation of geometry along the axis 
of the hole, there is no need for more than one subelement in that direction. The 
simplest conceivable representation of hole geometry is a single subelement with 
the volume integration for it's stiffness matrix extending only over the volume of 
the actual material present. With the eight corners as integration points. Figure 
20 shows the Young's modulus and shear modulus predictions as a function of hole 
volume fraction for an epoxy unit cell. Similar estimates may be obtained using 
3x3x2 and 4x4x2 integration grids depicted in Figure 21. However, the single 
eight-node subcell is not capable of yielding better stiffness estimates than the 
rule of mixtures (for elements in parallel) for this example, irrespective of 
integration grid refinements. 

Now consider a model of the same microstructure with more subcells of smaller 
size. The unit cell model in Figure 22 contains 12 subcells. However, only two 
are elastically different, the square subcells at the four corners and the eight 
rectangular ones between the corners. Stiffness estimates for this model are 
shown in Figure 23. Both 2x2x2 and 3x3x2 grids are used for subcell stiffness 
integrations. There is a small difference in results for the 2x2x2 and 3x3x2 
grids but there is no major improvement in accuracy with this or further grid 
refinement. However, a large improvement has resulted from the use of smaller, 
more numerous subelements (Compare Figures 20 and 23) . The standard of comparison 
in Figure 23 is the strength of materials calculation for stretching and bending 
of two orthogonal sets of parallel plates, welded together at their lines of 
intersection. This model is not an exact solution but is very accurate for hole 
volume fractions of more than 50%. 
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For high void fractions, moduli prediction of the fabric analysis and predictions 
of the strengths of materials model are generally in agreement; except for the 

shear modulus case. These predictions differ by almost a constant factor of 2. 

The reason for this difference is that the G shear modulus is dependent on 

xy 

bending stiffness of the thin cellular walls between the square holes. The eight 
node hexahedral element used to represent these cell walls is deficient in its 
ability to model linear bending strain distributions. 


33 



9. WOVEN FABRIC APPU:CA!r.lQNS 

This section considers application of the proceeding analysis to woven fabric 

microgeometry . First consider the most common weave construction, plain weave. 

Idealizing the geometry is traditionally the first step in analyzing any 

structure. For fabrics this usually implies replacing tow bundles with elastic 

tubes of constant cross section. In textile mechanics these interwoven tubes are 

referred to as channel models (Ref. 18) . In composite mechanics these tubes are 

assumed to consist of unidirectional orthotropic material with a plane of isotropy 

normal to the principal axis of each tube. Figure 24A contains one such 

idealization in which tow cross-sections are assumed to have a diamond shape. 

Dissimilar materials in the unit cell are separated by a series of flat planar 

areas. As mentioned previously, if warp and fill materials and their geometries 

are identical, and if a rule of four subcells per yarn crossover is adopted, then 

only one master subcell plus 16 transformations of that subcell are required to 

model the unit cell, (Figure 24B,C) . If a 2x2x2 integration grid is superimposed 

on the idealized subcell microgeometry, then the small volumes of dissimilar 

materials associated with each integration point can be visualized, characterized, 

and tabulated from geometric considerations (Figures 25 and Table 2). Table 2 
presents the two spherical fiber angles ((J)^,(J) 2 ) and eight material volume 

fractions (u , u , ‘ * * "O, f.) associated with each octant of volume surrounding each 

12 o 

integration point for each of the three constituent materials, namely, the 

unidirectional warp composite (1), the unidirectional fill composite (2), and the 
bulk matrix (3). The subscripts 1 through 8 on the material volume fractions 

refer to the following x,y, z octants respectively; +++, ++-, -++, -+-, - 

-+, and as shown in Figure 12. 

The volume fraction of impregnated tow in this composite is 50%. The other 50% is 


34 



interstitial bulk matrix. If elastic properties of the transversely isotopic tow 

material have the same values as the unidirectional material of the previous 

section (Equation 8.1) and the bulk ^^latrix properties are 

E = 0.5 msi, G *= 0.185 msi, V = 0.35, (9.1) 

then the analytical predictions of the fabric reinforced composite properties are: 
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If this same graphite/epoxy composite was unwoven unidirectional material in the 
form of a thick blended 0/90 laminate with 50% bulk matrix material between the 


plies, the laminate properties would be: 
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The weave microgeometry is thus responsible for an approximate 22% reduction in 
Young's moduli E^ and E^, but a much smaller reduction in the in-plane shear 


moduli . 

This sample problem serves only as an example of the application of fabric 
analysis to a reinforcing geometry resembling a plain weave microstructure. A 
true fabric microstructure has more complex tow paths and tow cross-sectional 
variations than this example. However, the analysis is not limited by any of 
these geometric complications, as the next example illustrates. 

9.1 Realistic Plain Weave Model 


Now consider the microgeometry of Figure 26 for a graphite/epoxy, plain weave, 
reinforced composite which is magnified approximately 70X. This laminate is an 8- 
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ply composite consisting of T-300 untwisted tow yarns at 18 ends and 18 picks per 
inch. The matrix designation is 5208; average thickness per ply is 0.011 inches; 
fiber volume fraction is 66.7%. Closer examination of Figure 26 reveals that warp 
and fill fibers have a maximum angle of approximately ± 6^ with the middle plane 
of the fabric as they undulate over and under one another. Using the concept of a 
somewhat irregular unit cell and subcell of Figure 9 with a 2x2x2 integration 
network (the corners of the subcell are the integration points ) , Table 3 contains 

all necessary remaining geometric input data for computation of the subcell 
stiffness matrix. As before, 1)^ designates the fraction of volume associated with 

the i th octant of volume (surrounding an integration point) that is occupied by 
one of the three constituent materials. Octants containing both warp and fill 
tows are considered to contain two different orthbtropic materials. Bulk matrix 
is the third material. 

Subcell cross-section sketches and 3-D yarn bundle models, based on 
photomicrographs, assist in establishing the 1)^ quantities. For simplicity, fiber 

angles are obtained from observed values at integration points rather than from 
averages over neighborhoods of the integration points. There is no attempt in 
this example to force reinforcing geometry to coincide with any idealized model. 

The fiber volume fraction within a tow bundle is taken to be 75% based on 
electronic scanning of similar composite photomicrographs using only those areas 
within the tow cross sections. The void content of the composite is assumed to be 
zero. The volume fraction of fiber in the tow (75%) , times the volume fraction of 
impregnated tow within the composite must equal the total fiber volume fraction of 
64%. 
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The volume fraction of impregnated tow is thus 85% of the total volume.* The 
volume fraction of the composite occupied by unreinforced matrix material must 
therefore be 15%. This is necessary to arrive at an overall composite fiber 
volume fraction of 64% which is the approximate average measured value for all of 
the test materials considered in this section (as determined by acid digestion) . 

There is some difficulty associated with the assignment of principal moduli to the 

impregnated tow material. Neither existing test data or micromechanics provides 

reliable estimates. Test data from unidirectional material seldom extends into 

the 75% fiber volume fraction range. Also, unidirectional test data does not 

include the inevitable degradation to tow properties that result from the weaving 

and related fabric forming processes. The latter objection also applies to 

micromechanics estimates of moduli. A semi-empirical application of the rule of 

mixtures (for elements in parallel) described subsequently seems to provide the 
best basis for establishing of the impregnated tows. 

In Ref. 19 twelve different unidirectional graphite composite materials were 

characterized experimentally. Their measured longitudinal moduli were lower in 

eleven out of twelve cases from the rule of mixtures prediction. The mixtures 
rule overestimated measured moduli (E^) by almost 10% based on an average of the 

twelve materials. This shortfall in the measured E^ can be attributed in part to 

fiber loss, misalignment, and breakage in the unidirectional prepregging and 
curing processes. The weaving process is considerably more damaging than the 
unidirectional prepregging process. 

*For analysis purposes 64% fber volume fraction was used rather than 66.7% because it was desired to have a common 
basis of conparison for eadi of the different weaves in this corrdatioa 
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Thus it seems reasonable to anticipate a greater reduction in due to weaving 
than unidirectional tow placement. Another 10% reduction in might be 

appropriate to account for weaving factors. Therefore, if the fiber volume 
fraction of a graphite/epoxy unidirectional material, with a 20 msi longitudinal 
modulus, were increased from 65% to 75% the rule of mixtures would predict about a 
15% longitudinal modulus increase. However, the weaving reduction factor would 
decrease this gain to only 5% and the resulting longitudinal modulus of the 
impregnated tow material within the weave would be about 21 msi. A truly reliable 
alternative to such an estimate would be an experimental study that measured 
impregnated tow modulus before and after weaving. This was beyond the scope of 
this program. 

The weaving process often includes beaming, sizing, weaving, scouring, drying and 
packaging. Each of these steps abrades, damages and misaligns the reinforcing 
fibers to some degree. The floor of any weaving room is a testament to the 
degradation and loss of reinforcing material. Also, T300 fibers showed more 
evidence of property reduction as a result of the unidirectional prepregging 
process than most of the graphite fibers studied in Ref. 19. 

Only the longitudinal modulus of the tow composite is assumed to be significantly 
degraded by weaving. Micromechanics considerations indicate a 15% increase in the 
transverse Young's modulus and shear moduli are appropriate to a fiber volume 
fraction increase from 65% to 75%. The principal Poisson's ratios should decrease 
a few percent as the fiber content increases. 

The principal moduli values for the impregnated tow composite are thus estimated 
to be (using the Rule of Mixtures for elements in series and parallel): 
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\ 

= 21.0 msi V^2 * ^*23 

E 2 *= E 3 - 1.7 msi V ^3 - 0.30 \ <9. 2) 

^12 0.8 msi G 23 =* 0.8 msi 

y 

Matrix properties are equivalent to those in the prior example (Equation 9.1). 
The fabric analysis predictions are given in Table 4 along with the experimental 
data. Experimental values and photomicrographs were provided by NASA Langley 
Research Center and were previously published in Reference 12. 

The comparison between analysis and experiment is generally good. The small 
differences in Young's moduli would indicate that the assumptions regarding tow 
property reduction resulting from the weaving processes were reasonable. The 
experimental difference between the moduli in the warp and fill direction could be 
accounted for in two ways. There are possibly some small variations in undulation 
angles between warp and fill tows that are not in evidence in the small area 
samples that were subject to microscopic examination. Also, there is ample reason 
to believe that property damage due to weaving is not evenly distributed between 
the warp and fill tows. 

There is a major discrepancy between the analytical and measured in-plane 
Poisson's ratio. The experimental value is suspect in this case because it is 
much greater than other analytical predictions, similar graphite /epoxy data, and 
0/90 cross ply analysis and data. 

The reduction in the principal in-plane moduli due to the weave microstructure is 
about 5% based on a cross-ply unidirectional laminate with the same fiber volume. 

Consider the possible loss of accuracy resulting from doubling the subcell x,y 


39 



side lengths while reducing the number of subcells in the unit cell from 16 to 4. 
One subcell then represents one warp/fill yarn crossover. This increases the 
maximum subcell side length ratio from 5/4 to 5/2, which is not excessive for 
finite element analysis. If a 3x3x2 integration grid is applied to the subcell 
then all the fiber angle and material volume fraction data associated with each 
integration point carries over unchanged from the 16 subcell model to the 4 
subcell model. 

Use of this larger subcell leads to the following plain weave fabric reinforced 
composite moduli predictions: 



Small Subcell (msi) 

Large Subcell (msi) 

E , E 

X y 

9.25 

9.22 

E 

z 

1.65 

1.65 

G , G 
yz xz 

0.721 

0.720 

G 

xy 

0.699 

0.744 



Small Subcell 

Large Subcell 

^yz'^xz 

0.329 

0.333 

V 

'^xy 

0.031 

0.028 


Results from the two models are almost identical. 

9.2 Other Weaves 

Using the same fiber and matrix, NASA has made composite laminates in four weave 
patterns: plain, 2/2 Oxford, five harness satin and eight harness satin. All 

weaves have 18 ends and 18 picks per inch. Figure 27 shows photomicrographs of 
the latter three weave geometries after lamination. Different weave patterns 
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yield slight variations in fiber undulation angles; The maximum angle that satin 
weave fibers make with the plane of the fabric is approximately one degree smaller 
than the maximum plain weave fiber angle. Oxford weave angles are smaller than 
plain weave angles but larger than satin weave angles. Fiber volume fractions of 
the four weaves also vary slightly. In the analysis a fiber volume fraction of 
64% was maintained for all weaves. 

Figure 28 shows an Oxford weave unit cell and three possible subdivisions. The 
first possibility (Figure 28C), with four subcells per yarn crossing/ requires 32 
transformations of two master subcells to model the unit cell. This model leads 
to a matrix inversion of dimension 93 square; probably larger than warranted. The 
second possibility (Figure 28D), with two subcells per yarn crossover, leads to a 
matrix inversion of order 45. Only 16 transformations of a single master subcell 
are required with this model. A third possible model is shown in Figure 28E. It 
is a coarser subdivision than either previous model. Each subcell represents one 
warp/fill tow crossover. All three models are shown to illustrate the point that 
many variations are possible with this type of analysis. The fabric reinforced 
composite moduli prediction, based on the medium subcell division model, using a 
2x3x2 integration grid, are given in Table 4. 

Comparison of the Oxford weave and plain weave moduli predictions shows a slightly 
greater Young's moduli for the Oxford weave. This reflects the differences in 
yarn crossover angles and crossover frequency. Fiber volume fraction for the 
Oxford weave was 64% for the analysis and 62% for the experiment. 

Consider the five harness satin weave geometry shown in the photomicrograph of 
Figure 27 and the sketch of Figure 29. If a subcell division of the unit cell is 
based on a rule of four subcells per yarn crossing then the largest matrix 
inversion is of order 297. Use of one subcell per crossover reduces this 
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dimension to 72. Thus, the suboell division shown in Figure 29C was adopted. Two 
different master subcell stiffness matrices are required for this model. A 3x3x2 
integration grid is used on each subcell. Matrix and unidirectional properties 
used are equivalent to those used in the Oxford and plain weave models. Analytical 
and experimental moduli are given in Table 4. 

Analysis again predicts the trend toward higher in-plane moduli with decreasing 
density of yarn undulations. The frequency of undulations has little effect on 
shear moduli or Poisson's ratios. Both analytical and experimental fabric volume 
fractions of the five-harness satin weave are 64%. 

The in-plane Young's moduli correlation is not as good as was obtained on the 
plain weave or Oxford weave. This raises a question concerning the use of a 
constant unidirectional tow composite property reduction factor to account for tow 
damage in weaving. It would appear from the correlation that the amount of tow 
damage is a function of the weave style. The lower than expected analytical 
Young's moduli for the satin weave, with relatively few warp/fill crossovers per 
unit of fabric area, indicates a possible lower level of tow damage than was 
evidenced in the plain or Oxford weave forming process with many more warp/fill 
crossovers. The beat up process following pick yarn insertion could be more of a 
localized damage phenomenon in the vicinity of warp/fill crossovers than an 
overall damage mechanism. 

l 

The eight harness satin weave is shown in the photomicrograph of Figure 27 and 
sketched in Figure 30. One unit cell is shown in Figure 30B. Using four subcells 
per yarn crossing leads to a reduced stiffness matrix of 765 square. One subcell 
per yarn crossing gives a stiffness matrix of order 189. In addition, three 
different master subcells are required. Thus, it is of interest to consider a 
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cruder single subceli that includes four yarn crossovers. Only two different 
master subcells require consideration with such a model. This subdivision of the 
unit cell is shown in Figure 30C. The unit cell reduced stiffness matrix is of 
order 45. Use of the previous constituent material properties with the 5x5x2 
integration grid of Figure 30D yields the results given in Table 4 . 

The lower than expected analytical Young's moduli in the plane of the fabric for 
the eight harness satin reflects the same trend as was observed for the five 
harness satin and again suggests less weave daimage with fewer warp/fill 
crossovers . 

In general the correlation between analysis and experiment was satisfactory for 
most engineering applications. A linearly varying tow property reduction factor 
based upon the number of warp/fill crossovers per inch of tow would have very much 
improved the correlation. However, this correction should be verified more 
thoroghly before adoption. 

In summary, a general rule of "four subcells per yarn crossover with crude 
integration schemes and networks" is adequate for modeling conventional woven 
fabric reinforcing geometries. Larger subcells can be used with little compromise 
in accuracy if corresponding refinements are made in the integration network. One 
subcell per ply in the thickness direction seems adequate under the same 
circumstances. 
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XSL BRAIDE D . FABRI C AegLXCAIIQNS 

A variety of industrial braided fabrics have application as composite reinforcing 
materials. Their microgeometries are often similar to woven fabrics. The 
simplest 2-D braids (Figure 31A) are analogous to skewed plain weaves. 

One important braid characteristic is braid angle, i.e., the average angle (in the 
plane of the fabric) that yarns (tows) make with the fabric output (machine take- 
up) direction. A braid angle of ± 45^ corresponds to an orthogonal plain weave 
(although some microgeometry differences may result from differences in tow 
handling) . For analytical purposes ± 45^ plain braid and plain weave 
microgeometry can be considered equivalent. Furthermore, when the braiding tows 
are identical and have identical spacing, it is possible to isolate a unit cell 
that is a rectangular prism (Figure 31B) . This particular unit cell can be 
subdivided into four subcells (Figure 31D) , each of which can be obtained from the 
other by various coordinate transformations. Thus, it is necessary to obtain a 
stiffness matrix for only one master subcell in order to model the unit cell. The 
stiffness matrix for the subcell is obtained by the same method used for weaves. 
First, a 3x3x2 integration grid is superimposed on the subcell volume. Fractional 
values of octant volumes (surrounding each integration point) containing the three 
constituent materials (two sets of impregnated braid tows and bulk matrix 
insterstices) are given in Table 5. Two spherical angles describing local fiber 
directions at each of the 18 integration points are also tabulated. This data 
corresponds to a braid angle of ± 45*^. The unidirectional composite properties 
are given in equation (9.2). Bulk matrix properties are given in equation (9.1). 

Moduli predictions for the 64% fiber volume fraction braid and for the 64% fiber 
volume fraction plain weave from Section 9. should be close to equivalent after 
the plain weave moduli are transformed into the coordinate system of the braid. 
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The predicted moduli are; 

Transformed 

± 45^ Braid Plain Weave 

E ,E (msi) 

X y 

E (msi) 
z 

G (msi) 
xy 

V 

xy 

Small differences in the descriptions of the microgeometry of the master subcells 
would account for the small moduli differences. The reason for the large 
difference in the Poisson's ratio of the two models is not clear. 

Figure 32 contains a plot of variation in braid moduli as a function of braid 
angle. Reinforced and unreinforced material volumes associated with integration 
points are assumed to remain constant as the braid angle varies. Actually/ some 
of the braid angles are not possible to achieve without "jamming” of the two braid 
tow systems. Distortions of local geometry begin as the braid angle approaches 
these limits. Braid moduli plots are qualitatively similar to the corresponding 
plots for symmetric angle-ply laminates made from unidirectional materials. 

Braids with different sets of braid tows may be analyzed by the same procedure, 
but more than one master subcell is needed to build the unit cell. Different tow 
spacings in the two sets of braided tows would present a greater modelling 
problem. 



A triaxial braid reinforced composite, as shown in Figure 33, consists of three 
sets of tows, intertwined together. Besides the pair of conventional braider tows 
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the 0° or longitudinal "stuffer" tows are introduced into the 2-D braid through a 
set of stationary tow carriers around which the braider carriers pass during the 
fabric formation. The 0° tows remain essentially straight in the final fabric and 
lie in the machine take up direction. The braid pattern must be modified to make 
room for the 0° tows. Instead of each braid tow intertwining over and under each 
braid tow that it crosses, each braid tow is allowed to float alternately over and 
under each pair of braid tows it crosses. This relates the braid to a 2/2 twill 
weave in the same way that the previous braid relates to a plain weave. The 
triaxial braid introduces an additional 0° tow size parameter of choice. For the 
particular triaxial braid to be modeled the tows were all made from AS-4 graphite 
fibers. The braid angle was ± 70°. The braider tows contained 6000 strands per 
tow. The 0° longitudinal tows contained 18,000 strands. Both tows were untwisted. 
The composite fiber volume fraction was 52% as measured by acid digestion. The 
matrix matrial was Shell 1895 epoxy. 

Figure 34 shows the smallest rectangular unit cell for the triaxial braid. It 
contains two 0° stuffer tows and portions of four tows from each of the two sets of 
±70° braider tows. The spacing of the stuffer tows establish the width of the unit 
cell. The stuffer tow spacing was 4.4 tows per inch making the unit cell (2/4.4) 
= 0.455 inches wide. The braid angle of ±70° established the unit cell height of 
(1/2) (0 . 455) / (tan 70°) - 0.166 inches. The thickness per ply was 0.0275 inches. 
The laminate was five plies thick. Figure 35 shows two photomicrographs of the 
actual material which was made by the Boeing Co. and tested at NASA Langley 
Research Center (Ref. 20) . The braider tows had an average crossover angle of ±9°. 
The unit cell was divided into eight subcells, as shown in Figure 34. The eight 
subcells reduce to two master subcells whose stiffness matrices were assembled 
using a 2x2x2 integration network. The fiber volume fraction within a tow bundle 
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was taken to be 75%. Thus, the volume of unidirectional tow composite within the 
unit cell was (100) (52/75) » 70%. The remaining 30% is bulk interstitial or 
unreinforced matrix material. The same elastic properties that were attributed to 
the bulk matrix and the unidirectional tow composite in the prior braid analysis 
were applied to the triaxial braid analysis. The percentages of fiber in the 
three different tow directions were almost equal. Since the braid angle was close 
to 60 ° it is expected that the resulting composite will have close to quasi- 
isotropic laminate properties for the same volume fraction of the same fiber. The 
principal triaxial composite moduli from both the fabric analysis and test are, 
with reference to the coordinate system of Figure 34: 



Analysis 

Test 

Ejj (msi) 

7.41 

7.03 

Ey (msi) 

6.32 

6.31 

E^ (msi) 

1.75 

N.A. 

Gjjy (msi) 

1.80 

N.A. 

Vxy 

0.214 

0.190 

Vyx 

0.169 

0.183 


The analysis and test are in good agreement, particularly as relates to measured 
strains in a tow diection. 
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The foregoing braided and woven fabric composite analysis proceedures also apply 
to composites with knitted reinforcement. Low fiber content and the absence of 
long fiber floats in the microgeometry render this form of reinforcement 

inefficient for highly-loaded structural applications. However, knitting is used 
in nonstructural applications because of its ability to conform to complex 
. surfaces prior to curing and because of its availability in tubular form. 

One of the most common knit patterns is the Jersey (plain) knit. Jersey 

microgeometry and a rectangular unit cell (Ref. 21) are shown in Figure 5. This 
pattern resembles some forms of chain link fencing which have almost the same unit 
cell. One possible subdivision of the unit cell into subcells is shown in Figure 
36. This subdivision is convenient because all four subcells can be formed from 
transformations of one master subcell. A 2x2x2 integration network serves to form 
the stiffness matrix for the master subcell. Spherical fiber direction angles and 
composite and bulk matrix volume fractions associated with each integration point 
are given in Table 6. The fiber volume fraction within the impregnated tow is 

assumed to be 75% and the volume fraction of unidirectional material within the 
composite is assumed to be 33.3%. Thus, the overall fiber volume fraction of the 
composite is (0.75 x 0.33 x 100) = 25% and interstitial or unreinforced matrix 
volume fraction is 66.7%. Unidirectional composite and matrix constituent 
properties are the same as the previous braided fabric constituent properties. 

The predicted knit fabric reinforced composite moduli are as follows: 
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The in-plane composite Young's moduli are lower than (0/90) laminates of 
unidirectional material with the same fiber content. The in-plane shear modulus 
and Poisson's ratio are higher for the knit. 

In braids, weaves, and knits prestress can distort microgeometry within the unit 
cell and cause large changes in moduli of the composite material. Prestretching 
of knits in the layup process is almost unavoidable. This phenomena is one reason 
for the absence of reliable test data on knits. Figure 37 shows the effects of 
prestretching or elongating the unit cell. As in the case of braids, tow 
compaction may prevent some of these microgeometries from being achievable. 
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12. CONCLUDING REMARKS 

This report describes a method for applying general 3-D finite element analysis to 
predict all of the linear elastic constants of fabric reinforced composite 
laminates. The analysis presumes the reinforcing microstructure can be condensed 
into a small repeating rectangular unit cell. The generic composite structure is 
presumed to consist of an infinite array of these unit cells joined together at 
matching surfaces. The analysis of a single unit cell is sufficient to predict 
all the global elastic properties of the composite. The approach avoids some of 
the difficulties associated with matching of finite element node point locations 
and element boundaries to internal material boundaries. Partially reinforced 
elements are used, along with the general energy formula for their stiffness 
matrix. This formula accounts for material property variations within the 
element. This approach should work if element dimensions are small with respect 
to the tow cross-sectional dimensions and interstitial matrix volume dimensions. 
This report shows that for many reinforcing geometries it is not essential (within 
the requirements of engineering accuracy) that the element size be that small. 
With most common reinforcing weaves and with the eight-node hexahedral element, a 
general rule for sizing of elements is one element per ply in the thickness 
direction and one to four elements per tow crossover. This rule applies to braids 
as well. Its application to knits is suggested but unverified. 

Since most general purpose finite element codes are awkward to use with the 
boundary conditions associated with moduli prediction, a special purpose computer 
code was written to facilitate application to these problems. This code is much 
quicker and easier to use than general finite element codes. Appendices I, II and 
III contain a listing, discussion of input data and sample problem for this 
Fortran code. 
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One advantage to this analysis is that it easily adapts to most forms of 
reinforcing microgeometry. The method does not require difficult logic or 
multiple option choices and enables the calculation of the all 3-D elastic 
contents. Hence, this method is useful for estimating material property inputs 
required of impact damage analysis codes. Another advantage of this analysis is 
that it remains within the scope of finite elements, which is familair to many 
engineers. No depth of knowledge in textile structures is required. 

Examples of application of this analysis to weaves, braids and knits are included 
to familiarize readers with fabric microgeometry and modeling of yarn 
constructions, and to impart some appreciation for the effects of varying these 
construction parameters. It is possible, with this code, to generate a catalogue 
of moduli predictions based On microgeometry and constituent property variations. 
However, the goal of this report is simply to show that reinforced composite 
microgeometry is amenable to routine structural analysis and that this analysis is 
useful in the search for impact damage resistant composite designs. 

The analytical predictions of fabric reinforced composite moduli have been 
compared to test data for several weave and braid reinforced composite laminates. 
The predictions are largely within the desired range of accuracy required of most 
materials engineering applications. 

If a laminate consists of several different plies Of different fabrics then it is 
necessary to apply the fabric reinforced composite analysis to each different 
fabric. This ply level analysis is followed by a conventional laminate analysis 
in order to obtain the elastic properties of the complete laminate. Laminates 
that are made from many plies of the same fabric reinforcement, oriented at 
different angles with respect to each other require a single fabric reinforced 
composite analysis followed by a conventional laminate analysis. Thin laminate 
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response may only be approximated using this analysis. Corrections may be 
necessary to account for free surface effects and possible single ply bending and 
stretching/bending coupling effects. 

Historically very few special elements have been developed which contain more than 
one material within the elements (except for laminated plate and shell elements) . 
The potential for new 3-D elements exists within the context of this analysis. 
However, this work is beyond the scope of the present effort. 
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0 
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TABLE 1: UNIT CELL CORNER DISPLACEMENTS FOR UNIT STRAIN CASES 
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- 
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0 
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- 
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7 

0 
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- 
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- 
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2 
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0 
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0 
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0 
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- 
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3 

3 

0 
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- 

- 
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3 

4 

0 
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- 
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- 

- 
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5 

0 
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0 

- 
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-■ 

3 

6 

0 
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- 

- 

- 

- 

- 

0 

- 
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3 

7 

0 
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- 

- 

- 

- 

- 
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TABLE 3: PLAIN WEAVE INPUT DATA FROM PHOTOMICROGRAPHS 
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- 
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TABLE 5A: BRAID INPUT DATA (MATERIAL 
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TABLE 5B: BRAID INPUT DATA {MATERIAL 2) 
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FIGURE 1 : ROLE OF FABRIC MECHANICS IN COMPOSITE MECHANICS 
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Unit Cell 
(A) 


Design Plan 
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Unit Cell 
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FIGURE 3: PLAIN WEAVE DESIGN PLAN AND UNIT CELLS 



FIGURE 4: TWILL WEAVE DESIGN PLAN AND UNIT CELL 
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Jersey Knit Tricot Knit 


FIGURE 5: UNIT CELLS FOR WEAVES/ BRAIDS, AND KNITS 
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stack Bond Brick Plain Weave Stacking 


FIGURE 6: UNIT CELL ARRAY 
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FIGURE 8: IMAGE POINTS ON UPPER AND LOWER SURFACE OF UNIT CELL 
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Unit Cell 


FIGURE 9: PLAIN WEAVE GEOMETRY 
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Global Coordinates of Integration 


FIGURE 11: FIBER DIRECTION IN SPHERICAL COORDINATES ( ^ > <l>2 ^ 
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(B) Octants of Volume at an Integration Point 
FIGURE 12: INTEGRATION GRID GEOMETRY DETAILS 
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FIGURE 15; FLOW CHART FOR THE COMPUTER CODE 
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Subcell I 


FIGURE 16: 2-D RICKRACK COMPOSITE GEOMETRY 










FIGURE 17: 3-D RICKRACK COMPOSITE GEOMETRY 
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Hole Volume Fraction 


FIGURE 20: MODULI OF VOIDED MATERIAL BASED ON A SINGLE SUBCELL 
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FIGURE 21: REFINING THE INTEGRATION GRID 



(B) Subcells 



I II 


FIGURE 22: TWELVE SUBCELL MODEL 
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FIGURE 25: EXPLODED VIEW OF MASTER SUBCELL FROM FIG. 24(C) 
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(B) Filling Direction 


FIGURE 26: LAMINATED PLAIN WEAVE PHOTOMICROGRAPHS (70X) 
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5 Harness Satin Weave - 8 plies 



8 Harness Satin Weave - 8 plies 



FIGURE 27: LAMINATED OXFORD AND SATIN WEAVE PHOTOMICROGRAPHS 
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FIGURE 28: 2/2 OXFORD WEAVE MODELS 
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(C) Subcell Division 


FIGURE 29: FIVE HARNESS SATIN WEAVE MODEL 
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(D) Subcell Intergratlon Network 


FIGURE 30: EIGHT HARNESS SATIN WEAVE MODEL 
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FIGURE 32: VARIATION IN STIFFNESS WITH BRAID ANGLE 
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FIGURE 33: 2-D TRIAXIAL BRAID 
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FIGURE 34; 2-D TRIAXIAL BRAID MODEL 
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Braider Tow 


Stuffer Tow 



Transverse Cross-Section (20x) 


FIGURE 35: PHOTOMICROGRAPHS OF TRIAXIAL BRAID 
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(C) Subcell Assembly 


FIGURE 36: JERSEY (PLAIN) KNIT MICROGEOMETRY MODEL 
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FIGURE 37: EFFECT OF PRESTRAIN ON JERSEY KNIT COMPOSITE STIFFNESS 
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APPENDIX I 




Program Listing-Fortran 5 


PROGRAM FABNEW ( INPUT, OUTPUT , TAPE5=INPUT, TAPE6=OUTPUT) 
C 

C cccccccccccccccccccccccccccccccccccccccccccccccccc 

c cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C CC PROGRAM TO COMPUTE THE 3-D ELASTIC MODULI OF CC 

C CC A UNIT CELL OF FIBER REINFORCED COMPOSITE CC 

C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C SPECIFICATION STATEMENTS 

C 

REAL MU,MIXYZ,MIYYZ,MIZY2,MIXXZ,MIYXZ,MIZX2 
REAL MIXXY,MIYXY,MI2XY 

REAL MP(6,6) ,KS(24,24) ,KMOD(24,24) ,SK(24,24) 

DIMENSION KA(7) , LC(24) , VF(8) 

DIMENSION PR(6,6) ,SS(6,6) ,T(6,6) ,BIG(6,6) 

DIMENSION BM(6,24) ,DB(6,24) 

DIMENSION GM(6,9) ,DG(6,9) ,GDG(9,9) 

DIMENSION GB(9,24) ,BDG(24,9) 

PARAMETER (MM=10 , MMM=100 , NNN==2 16 ) 

DIMENSION PROP (MM, 6) ,DX(MM) , DY (MM) , DZ (MM) 

DIMENSION NBL(MM,MM) , FDX(MM+1) ,FDY(MM+1) 

DIMENSION FB(NNN,7) ,FS(MMM,7) , FT ( 6 , 6) , TF ( 6 , 6) 

DIMENSION IP(MMM) ,WK(MMM) ,UVW(NNN,6) 

REAL KB(NNN,NNN) ,KM(MMM,NNN) ,KN(MMM,MMM) 

INTEGER QT,QTR(MM,MM) 

LOGICAL LX,LY 

BUILT IN MATERIAL PROPERTY DATA 

DATA (MP(1,I) ,I=1,6)/30.E6,1.5E6, .25, .35, .7E6, ,7E6/, 

1 (MP(2,I) , 1=1, 6)/2 5.E6, 1. 5E6, .25, .35, .7E6, .7E6/, 

2 (MP(3,I) ,I=1,6)/20.E6,1.5E6, .3, .45, .7E6, .7E6/, 

3 (MP(4,I) ,1=1, 6)/10.E6, 1.E6, .3, . 4 , 2 . E5 , 2 . E5/ , 

4 (MP(5,I) ,1=1,6)/ l.E7,l.E7, .25, .25,4.E6,4.E6/ 

5 (MP(6,I) ,1=1,6)/ .5E6, .5E6, .35, .35, .185E6, .185E6/ 
DATA (KA(I) ,I=1,7)/10,6,6,1,6,0,0/ 

INITIALIZE VARIABLES 

ISYM=0 

DO 10 1=1, MMM 
DO 5 J=l,7 
5 FS(I,J)=0.0 
DO 10 J=1,MMM 
10 KN(I,J)=0.0 
DO 15 1=1,6 
DO 15 J=l,6 
TF(I,J)=0.0 
15 FT(I,J)=0.0 
WRITE (6, 9100) 

READ(5,9030) NM 

MATERIAL PROPERTY DATA INPUT 

DO 18 1=1, NM 
WRITE(6,9180) 
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READ(5,9030) M 
IF(M,GT.6) THEN 
WRITE(6, 9120) 

READ(5,9010) PROP(I,l) 
WRITE(6,9130) 

READ(5,9010) PROP(I,2) 

WRITE(6, 9140) 

READ (5, 9010) PROP (I, 3) 

WRITE (6, 9 150) 

READ(5,9010) PROP(I,4) 

WRITE (6, 9160) 

READ{5, 9010) PROP(I,5) 

WRITE(6, 9170) 

READ(5,9010) PROP(I,6) 

END IF 

IF(M.LE.6) THEN 
DO 17 J=l,6 

17 PROP(I, J)=MP(M, J) 

END IF 

18 CONTINUE 
WRITE(6, 9190) 

DO 19 I“1,NM 

19 WRITE(6,9020) ( PROP ( I , J) , J=1 , 6 ) 

WRITE(6,9110) 

READ(5,9030) NB 
WRITE (6, 9080) 

READ(5,9030) NBX 
WRITE (6, 9090) 

READ (5, 9030) NBY 

NP= ( (NBY+1) ★ (NBX+1) ) *6 

IJ«3* (NBX*NBY~1) 

WRITE(6, 9470) 

READ(5,9000) C 
WRITE(6, 9440) 

READ (5, 9000) XL 
FDX(1)=0.0 
FDX(NBX+1) =100.0 
IF(NBX.LE. 1) GO TO 25 
DO 22 I=1,NBX-1 
WRITE(6, 9460) I+l 
22 READ(5,9000) FDX(I+1) 

25 WRITE(6, 9450) 

READ (5, 9000) YL 
FDY(1)=0.0 
FDY(NBY+1)=100.0 
IF(NBY. LE. 1) GO TO 35 
DO 30 I=1/NBY-1 
WRITE(6, 9460) 1+1 
30 READ(5,9000) FDY(I+1) 

35 DO 50 1=1, NP 
DO 40 J=l,7 
40 FB(I, J)=0.0 
DO 50 J=1,NP 
50 KB(I,J)=0.0 
DO 70 1=1, NBX 
DO 70 J=1,NBY 
WRITE(6, 9060) I,J 
READ(5,9030) NBL(I,J) 
WRITE(6,9070) I,J 
READ(5,9030) QTR(I,J) 
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70 CONTINUE 

WRITE (6, 9050) 

READ(5,9040) LX 
IF (LX) IJ=IJ+(NBY+1) 

WRITE(6,9055) 

READ (5, 9040) LY 
IF(LY) IJ=IJ+(NBX+1) 

I COUNTED 
350 CONTINUE 
C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C OUTER DO LOOP ON NO. OF SUBCELLS IN UNIT CELL BEGINS C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INPUT BLOCK GRID GEOMETRY 

WRITE(6,9300) 

WRITE (6, 9 2 00) ICOUNT+1 
WRITE (6, 93 00) 

WRITE(6, 9210) 

READ(5,9000) A 
AA=0.5*A 
WRITE(6,9240) 

READ(5,9030) NX 
DX(1)=0. 0 
DX(NX)=A 

IF (NX.LE.2) GOTO 460 
DO 450 1=2, NX-1 
WRITE (6, 9270) I 
READ(5,9010) DX(I) 

450 DX(I)=DX(I) *A/100.0 
460 CONTINUE 

WRITE(6,9220) 

READ (5, 9000) B 
BB=0.5*B 
WRITE(6,9250) 

READ(5,9030) NY 
DY(1)=0.0 
DY(NY)=B 

IF (NY.LE.2) GOTO 510 
DO 500 1=2 , NY-1 
WRITE(6,9280) I 
READ(5,9010) DY(I) 

500 DY(I)«DY(I) *B/100,0 
510 CONTINUE 

WRITE (6, 9230) 

READ (5, 9000) HT 
CC=0.5*C 

IF(.NOT.(HT.EQ.C)) WRITE(6,9330) 

WRITE (6, 92 60) 

READ(5,9030) NZ 
DZ(1)=0,0 
DZ(NZ)=C 

IF (NZ. LE.2) GOTO 560 
DO 550 1=2 ,NZ-1 
WRITE(6, 9290) I 
READ(5,9010) DZ(I) 

550 DZ(I)=DZ(I)*C/100.0 
560 CONTINUE 
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GENERAL MATERIALS INFORMATION FOR THE SUBCELL 

WRITE (6, 9310) ICOUNT+1 
READ (5, 9030) NMAT 
MCOUNT«0 
DO 570 1=1,24 
DO 570 J=l,24 
570 KS(I,J)*0.0 

DO 580 1=1,24 
DO 580 J=l,9 
580 BDG(I,J)=0.0 
DO 590 1=1,9 
DO 590 J=l,9 
590 GDG(I,J)=0.0 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
INNER DO LOOP ON MATLS . WITHIN SUBCELL BEGINS HERE C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
600 WRITE(6,9320) 

READ(5,9030) MN 

INPUT 8 MATL. VOL. FRACTIONS AT EACH INTEGRATION POINT 


DO 2000 1=1, NX 
DO 2000 J=1,NY 
DO 2000 K=1,NZ 
DO 800 LL=1,8 
800 VF(LL)«0.0 

WRITE (6, 9480) 

READ(5,9000) A1 
WRITE(6, 9490) 

READ (5, 9000) A2 
850 WRITE(6, 9420) I,J,K 
READ(5,9030) L 

IF( (L.LT. 0) .OR. (L.GT.8) ) GO TO 860 
WRITE(6,9430) L 
READ(5,9000) VF(L) 

GO TO 850 
860 CONTINUE 

IF(I.EQ.NX) THEN 

XP=0 . 0 

ELSE 

XP=(DX(I+1)-DX(I) )/2.0 
END IF 

IF(I.EQ.l) THEN 

XM=0. 0 

ELSE 

XM=(DX(I) -DX(I-l) )/2.0 
END IF 

IF(J.EQ.NY) THEN 

YP=0.0 

ELSE 

YP=(DY(J+1)-DY(J) )/2.0 
END IF 

IF(J.EQ.l) THEN 
YM=0. 0 

YM=(DY(J)-DY(J-1) )/2.0 
END IF 

IF(K.EQ.NZ) THEN 
ZP=0.0 
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ELSE 

ZP=(DZ(K+1)-DZ(K) )/2.0 
END IF 

IF(K.EQ.l) THEN 

ZM=0.0 

ELSE 

ZH=‘ (DZ (K) -DZ (K-1) ) /2 . 0 
END IF 

V1=VF(1) *XP*YP*ZP 
V2=VF(2) *XP*YP*ZM 
V3=VF(3) *XP*YM*ZP 
V4=VF(4) *XP*YM*ZM 
V5~VF(5) *XM*YP*ZP 
V6=VF(6) *XM*YP*ZM 
V7=VF (7) *XM*YM*ZP 
V8=VF ( 8 ) *XM*YM*ZM 
VOL=Vl+V2+V3+V4+V5+V6+V7+V8 

GET THE STRAIN TRANSFORMATION MATRIX (T) 

S1«=SIND(A1) 

S2=SIND{A2) 

Cl=COSD(Al) 

C2=COSD(A2) 

S1S=S1*S1 

S2S=S2*S2 

C1S=C1*C1 

C2S=C2*C2 

SC1=S1*C1 

SC2=S2*C2 

T(l, 1)=C1S*C2S 

T(1,2)=S1S*C2S 

T(l, 3)=S2S 

T(1,4)=SC1*C2S 

T(l, 5)=S1*SC2 

T(1,6)=C1*SC2 

T(2,1)=S1S 

T(2,2)=C1S 

T(2,3)=0.0 

T(2,4)=“SC1 

T(2,5)=0.0 

T(2,6)=0. 0 

T(3,1)=C1S*S2S 

T(3,2)=S1S*S2S 

T(3 , 3)=C2S 

T(3,4)=SC1*S2S 

T(3,5)=-S1*SC2 

T(3, 6)=-Cl*SC2 

T(4, 1)=~2 .0*SC1*C2 

T(4,2)=~T(4,1) 

T(4,3)=0.0 
T(4,4)=(C1S-S1S) *C2 
T(4,5)=C1*S2 
T(4,6)=-S1*S2 
T(5, 1)=2 . 0*SC1*S2 
T(5,2)— T(5,l) 

T(5,3)=0.0 
T(5,4)=-(C1S-S1S) *S2 
T(5,5)=C1*C2 
T(5,6)=-S1*C2 
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T(6, l)=-2 . 0*C1S*SC2 
T(6 , 2) «^2 . 0*S1S*SC2 
T(6,3)=2.0*SC2 
T ( 6 , 4 ) «-SCl*SC2*2 . 0 
T(6,5)=S1*(C2S-S2S) 

T(6,6)»C1*(C2S-S2S) 

GET STRESS-STRAIN MATRIX (SS) IN MATL. COORD. SYSTEM 

DO 1100 11=1,6 
DO 1100 JJ=1,6 
100 SS (II, JJ)=0.0 
MG=MN 

MU=PROP(MC, 4) 

R-PROP(MC, 3) *PROP(MC,3) *PR0P(MC, 2)/PR0P(MC, 1) 

D= ( 1 . 0+MU) * ( 1 . O-MU-2 . 0*R) 

IF(D.LE.O.O) D=1.0 

SS ( 1 , 1 ) =PROP (MC , 1 ) * ( 1 . 0-MU*MU) /D 
SS(l,2)=PROP(MC,2) *PROP(MC,3) *(1.0+MU)/D 
SS(1,3)=SS(1,2) 

SS(2,1)=SS(1,2) 

SS(2,2)=PROP(MC,2)*(1.0-R)/D 
SS(2,3)=PROP(MC,2) *(MU+R)/D 
SS(3,1)=SS(1,3) 

SS(3,2)=SS(2,3) 

SS(3,3)=SS(2,2) 

SS(4,4)=PROP(MC,5) 

SS(5,5)=PROP(MC,6) 

SS (6,6)=PROP(MC,5) 


GET STRESS-STRAIN MATRIX (SS) IN COORDS. OF THE SUBCELL 

DO 1200 11=1,6 
DO 1200 JJ=1, 6 
SUM=0.0 

DO 1150 KK=1,6 

1150 SUM=SUM+SS(II,KK) *T(KK,JJ) 

1200 PR(II, JJ)=SUM 
DO 1225 11=1,3 
DO 1225 JJ=4,6 
T(II,JJ)=2.0*T(II, JJ) 

1225 T(JJ, II)=0.5*T(JJ, II) 

CALL MATOPS ( KA , T , DET , DDDD ) 

DO 1300 11=1,6 
DO 1300 JJ=1,6 
SUM=0 .0 

DO 1250 KK=1,6 

1250 SUM=SUM+T(II,KK)*PR(KK,JJ) 

1300 BIG(II,JJ)=SUM*VOL 
1350 CONTINUE 

SUBCELL GEOMETRY CALCULATIONS 

DO 1400 11=1,6 
DO 1400 JJ=1,9 
1400 GM(II, JJ)=0.0 
DO 1500 11=1,6 
DO 1500 JJ=1,24 
1500 BM(II, JJ)=0.0 
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FORM 8 NODE FINITE ELEMENT STIFFNESS MATRIX FOR SUBCELL 


X=DX(I)-AA 

Y=DY(J)-BB 

Z«DZ(K)~CC 

X1=»“(0. 125/AA) * (1.0-Y/BB) * (1. 0-Z/CC) 
Y1=-(0.125/BB) *(1.0“X/AA) *(1.0-Z/CC) 
Z1=-(0.125/CC) * (1.0-X/AA) * (1.0-Y/BB) 
X2="(0.125/AA) *(1,0~Y/BB) *(1.0+Z/CC) 
Y2=-(0. 125/BB) *(1.0-X/AA) *(1.0+Z/CC) 
Z2=+(0. 125/CC) *( 1.0-X/AA) ★( 1.0-Y/BB) 
X3=-(0.125/AA) *(1.0+Y/BB) *(1. 0-Z/CC) 
Y3=+(0.125/BB) *( 1.0-X/AA) *(1. 0-Z/CC) 
Z3=*-(0.125/CC) *(1. 0-X/AA) * (1.0+Y/BB) 
X4=-(0. 125/AA) *(1.0+Y/BB) * (1.0+Z/CC) 
Y4=+ ( 0 . 125/BB) *( 1 . 0-X/AA) *( 1 . 0+Z/CC) 
Z4=+(0. 125/CC) ★ (1. 0-X/AA) * (1. 0+Y/BB) 
X5=+ (0 . 125/AA) * (1. 0-Y/BB) * ( 1 . 0-Z/CC) 
Y5=-(0.125/BB) *(1.0+X/AA) * (1. 0-Z/CC) 
Z5=»-(0. 125/CC) *(1.0+X/AA) * (1.0-Y/BB) 
X6=+(0. 125/AA) *( 1.0-Y/BB) * (1.0+Z/CC) 
Y6=-(0.125/BB) *(1.0+X/AA) * (1.0+Z/CC) 
Z6=+ (0 . 125/CC) * ( 1 . 0+X/AA) * ( 1 . 0-Y/BB) 
X7=+ (0 . 125/AA) * ( 1 . 0+Y/BB) * ( 1 . 0-Z/CC) 
Y7=+(0.125/BB) *(1. 0+X/AA) * (1. 0-Z/CC) 
Z7=- (0. 125/CC) *(1. 0+X/AA) * (1.0+Y/BB) 
X8=+(0. 125/AA) *(1.0+Y/BB) * (1.0+Z/CC) 
Y8='+(0.125/BB) *(1. 0+X/AA) * (1.0+Z/CC) 
Z8=+(0. 125/CC) + (1. 0+X/AA) ★ (1.0+Y/BB) 
XX=~2.0*X/(AA*AA) 

YY=-2.0*Y/(BB*BB) 

ZZ=-2.0*Z/(CC*CC) 

GM(1, 1)-XX 

GM(2,5)=YY 

GM(3,9)=ZZ 

GM(4,2)=XX 

GM(4,4)=YY 

GM ( 5 , 6 ) =YY 

GM(5,8)=ZZ 

GM(6,3)=XX 

GM(6,7)=ZZ 

BM(1,1)=X1 

BM(2,2)=Y1 

BM(3 ,3)=Z1 

BM(1,4)=X2 

BM(2,5)=Y2 

BM(3 ,6)=Z2 

BM(1,7)=X3 

BM(2,8)=Y3 

BM(3,9)=Z3 

BM(1,10)=X4 

BM(2,11)=Y4 

BM(3,12)=Z4 

BM(1, 13)=X5 

BM(2,14)=>Y5 


BM(3, 15)“Z5 


BM(1, 16)*X6 


BM(2,17)=Y6 


BM(3, 18)=Z6 
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BM(1,19)-X7 

BM(2,20)=Y7 

BM(3,21)=Z7 

BM(1»22)»X8 

BM(2,23)=*Y8 

BM(3,24)-Z8 

BM(4,1)=Y1 

BM(5,2)=Z1 

BM(6, 3)=X1 

BM(4 ,4)=Y2 

BM(5, 5)“Z2 

BM(6,6)»X2 

BM(4,7)=Y3 

BM(5,8)=Z3 

BM(6,9)=X3 

BM(4, 10)=Y4 

BM(5,11)=Z4 

BM(6,12)~X4 

BM(4,13)=*Y5 

BM(5,14)=Z5 

BM(6,15)»X5 

BM(4 , 16)-Y6 

BM(5,17)^Z6 

BM(6,18)=X6 

BM(4,19)-Y7 

BM(5,20)=Z7 

BM(6,21)=X7 

BM(4,22)*=Y8 

BM(5, 23)=Z8 

BM(6,24)==X8 

BM(6,1)=Z1 

BM(4,2)=X1 

BM(5,3)=Y1 

BM(6,4)=Z2 

BM(4,5)=X2 

BM(5, 6)=Y2 

BM(6,7)=Z3 

BM ( 4 , 8 ) =X3 

BM(5,9)=»Y3 

BM(6,10)=Z4 

BM(4 , 11)=X4 

BM(5,12)=Y4 

BM(6,13)=Z5 

BM(4 , 14 ) =X5 

BM(5,15)=*Y5 

BM(6,16)=Z6 

BM(4,17)=X6 

BM(5,18)=Y6 

BM(6,19)=Z7 

BM(4,20)=X7 

BM(5,21)=Y7 

BM(6,.22)=Z8 

BM(4,23)=X8 

BM(5,24)=Y8 

DO 1600 11=1,6 

DO 1600 JJ=1,24 

DB(II, JJ)=0. 0 

DO 1600 KK=1,6 

1600 DB(II,JJ)=DB(II,JJ)+BIG(II,KK)*BM(KK,JJ) 
DO 1700 11=1,6 
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DG(II,1)=XX*BIG(II,1) 

DG ( II , 2 ) =XX*BIG ( II , 4 ) 

DG(II,3)=XX*BIG(II,6) 

DG(II,4)»YY*BIG(II,4) 

DG(II,5)=YY*BIG(II,2) 

DG(II,6)=YY*BIG(II,5) 

DG(II,7)»Z2*BIG(II,6) 

DG(II,8)=ZZ*BIG(II,5) 

1700 DG(II,9)=*ZZ*BIG(II,3) 

DO 1750 11=1,9 
DO 1750 JJ=1,9 
DO 1750 KK=1, 6 

1750 GDG(II, JJ)=GDG(II,JJ)+GM(KK,II) *DG(KK,JJ) 

DO 1800 11=1,24 
DO 1800 JJ=1,9 
DO 1800 KK=1,6 

1800 BDG(II, JJ)=BDG(II, JJ)+BM(KK,II) *DG(KK,JJ) 

DO 1900 11=1,24 
DO 1900 JJ=1,24 
DO 1900 KK=1,6 

1900 KS (II, JJ)=KS(II, JJ)+BM(KK, II) *DB(KK, JJ) 

2000 CONTINUE 

MCOUNT=MCOUNT+ 1 

IF(MCOUNT, LT.NMAT) GO TO 600 

CALL INV(GDG) 

DO 2020 1=1,9 
DO 2020 J=l,24 
GB(I,J)=0.0 
DO 2020 K=l,9 

2020 GB(I, J)=GB(I, J)+GDG(I,K) *BDG(J,K) 

DO 2030 1=1,24 
DO 2030 J=l,24 
KMOD(I, J) =0.0 
DO 2030 K=l,9 

2 03 0 KMOD(I, J)=KMOD(I, J)+BDG(I,K) *GB(K, J) 

DO 2040 1=1,24 
DO 2040 J=l,24 

2040 KS(I,J)=KS(I,J)-KMOD(I,J) 

C 

C PUT SMALL STIF. MATRIX (KS) INTO BIG STIF. MATRIX (KB) 
C 

ICOUNT=ICOUNT+l 
DO 2400 1=1, NBX 
DO 2400 J=1,NBY 
QT=QTR(I,J) 

NK=NBL(I,J) 

IF(.NOT.(NK.EQ.ICOUNT)) GO TO 2400 
IF(QT.LT.O) GO TO 2400 
DO 2050 11=1,24 
DO 2050 JJ=1,24 

2050 SK(II, JJ)=KS(II, JJ) 

IF(QT.EQ.O) GO TO 2100 
CALL TRNS(QT,SK) 

2100 LC(1)=( (NBY+1) * (I-l)+J-l) *6+1 
LC(13)=LC(1)+(NBY+1)*6 
DO 2200 K=2,12 
LC(K)=LC(K-1)+1 

2200 LC(K+12)=LC(K+11)+1 
DO 2300 11=1,24 
DO 2300 JJ=1, 24 
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III=LC(II) 

JJJ*LC(JJ) 

2300 JJJ) =KB(III, JJJ)+SK(II, JJ) 

2400 CONTINUE 

IF(ICOUNT.LT.NB) GO TO 350 
GO TO 2460 
2440 ISYM=1 

DO 2445 1=1, IJ 
DO 2445 J=1,IJ 
2445 KN(I,J)=0.0 

DO 2447 1=1, IJ 
2447 KN(I,I)=1. 0 

IF(LX) IJ=IJ-2*NBY+1 
IF(LY) IJ=IJ-2*NBX+1 
DO 2450 1=1, NP 
DO 2450 J=l,6 
2450 FB(I,J)=0.0 
2460 CONTINUE 

CALC. DISP. VECTORS FOR 6 HOMOGENEOUS UNIT STRAIN CASES 


DO 2500 I=1,NBX+1 

DO 2500 J=1,NBY+1 

K=( (I-l) *(NBY+1)+(J~1) ) *6+1 

FB(K, 1)=FDX(I) *XL/100.0 

FB(K+3,1)=FB(K,1) 

FB(K+1,2)=FDY(J) *YL/100.0 
FB(K+4,2)=FB(K+1,2) 

FB(K+5,3)=C 

FB(K,4)=FB(K+1,2) 

FB(K+3,4)=FB(K,4) 

FB(K+4,5)=C 
:500 FB(K+3,6)=C 

DO 2525 1=1, NP 
DO 2525 J=l,6 
525 UVW(I, J)=FB(I, J) 

USE FORCE B.C.S TO ELIMINATE INTERNAL SURFACE FORCES 
IN=0 

IF( (NBX.LE.l) .OR. (NBY.LE. 1) )GO TO 2950 

DO 2900 1=2, NBX 

DO 2900 J=2,NBY 

L=( (I-l) *(NBY+1)+(J-1) ) *6 

DO 2850 K=1,NP 

KM(IN+1,K)=KB(L+1,K)+KB(L+4,K) 

KM(IN+2,K)=KB(L+2 ,K)+KB(L+5,K) 

850 KM(IN+3 ,K)=KB(L+3 ,K)+KB(L+6,K) 

IN=IN+3 
900 CONTINUE 

USE FORCE B.C.S TO ELIMINATE FORCE AT SIDE NODES 

2950 IF(NBX.LE.l) GO TO 3150 
DO 3100 1=2, NBX 
L=6*(I~1) *(NBY+1) 

LL=L+NBY*6 
IF(LY) GO TO 3050 
DO 3000 K=1,NP 

KM ( IN+1 , K) =KB ( L+1 , K) +KB ( L+4 , K) +KB ( LL+1 , K) +KB ( LL+4 , K) 
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KM ( IN+2 , K) =KB ( L+2 , K) +KB < L+5 , K) +KB ( LL+2 , K) +KB ( LL+5 , K) 
3000 KM ( IN+3 , K) =KB ( L+3 , K) +KB ( L+6 , K) +KB (LL+3 , K) +KB ( LL+6 , K) 
IN=IN+3 
GO TO 3100 
3050 CONTINUE 

IF(ISYM.EQ.l) GO TO 3080 
DO 3075 K=1,NP 

KM ( I N+ 1 , K ) = K B { L+ 1 , K ) + K B ( L+ 4 , K ) 
KM(IN+2,K)=KB(L+3,K)+KB(L+6,K) 

KM ( IN+3 , K) =KB ( LL+1 , K) +KB ( LL+4 , K) 

3075 KM(IN+4,K)=KB(LL+3,K)+KB(LL+6,K) 

IN=IN+4 
GO TO 3090 
3080 CONTINUE 

DO 3085 K=1,NP 

KM(IN+1,K)“KB(L+2,K)+KB(L+5,K) 

3085 KM(IN+2,K)=KB(LL+2,K)+KB(LL+5,K) 

IN»IN+2 
3090 CONTINUE 
3100 CONTINUE 

3150 IF(NBY.LE.l) GO TO 3305 
DO 3300 J=2 ,NBY 
L*6*(J-1) 

LL=L+NBX* (NBY+1) *6 
IF (LX) GO TO 3250 
DO 3200 K=1,NP 

KM ( IN+1 , K) =KB ( L+ 1 , K) +KB ( L+4 , K) +KB (LL+1 , K) +KB ( LL+4 , K) 
KM (IN+2 , K) =KB (L+2 , K) +KB(L+5 , K) +KB (LL+2 , K) +KB(LL+5 , K) 
3200 KM(IN+3 ,K)=KB( L+3, K)+KB( L+6, K)+KB( LL+3, K)+KB(LL+6,K) 
IN=IN+3 
GO TO 3300 
3250 CONTINUE 

IF(ISYM.EQ. 1) GO TO 3280 
DO 3275 K=1,NP 

KM (IN+1 , K) =KB (L+2 , K) +KB (L+5 , K) 

KM ( IN+2 , K) =KB ( L+3 , K) +KB ( L+6 , K) 

KM ( IN+3 , K) =KB (LL+2 , K) +KB (LL+5 , K) 

3275 KM(IN+4,K)=KB(LL+3,K)+KB(LL+6,K) 

IN=IN+4 
GO TO 3290 
3280 CONTINUE 

DO 3285 K*=1,NP 

KM(IN+l,K)=KB(L+l,K)+KB(L+4 ,K) 

3285 KM ( IN+2 , K) =KB ( LL+1 , K) +KB ( LL+4 , K) 

IN=IN+2 
3290 CONTINUE 
3300 CONTINUE 
3305 CONTINUE 

IF(.NOT.LY) GO TO 3315 
L=6*NBY 

LL=L+6*NBX* (NBY+1) 

IF(ISYM.EQ. 1) GO TO 3311 
DO 3310 K=1,NP 

KM (IN+1, K)=KB(L+1,K)+KB( L+4, K)+KB( LL+1, K)+KB (LL+4 ,K) 

3310 KM ( IN+2 , K) =KB (L+3 , K) +KB ( L+6 , K) +KB (LL+3 , K) +KB (LL+6 , K) 
IN=IN+2 

GO TO 3314 

3311 CONTINUE 

DO 3312 K=1,NP 

3312 KM ( IN+1 , K) =KB ( L+2 , K) +KB ( L+5 , K) +KB ( LL+2 , K) +KB ( LL+5 , K) 
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IN=IN+1 

3314 CONTINUE 

3315 CONTINUE 
IF(.NOT.LX) GO TO 3325 
L=6*NBX* (NBY+1) 

LL=L+6*NBY 

IF(ISYM.EQ.l) GO TO 3321 
DO 3320 K=1,NP 

KM ( IN+1 , K) =KB ( L+2 , K) +KB ( L+5 , K) +KB ( LL+2 , K) +KB ( LL+5 , K) 

3320 KM(IN+2,K)=KB(L+3,K)+KB(L+6,K)+KB(LL+3,K)+KB(LL+6,K) 
IN=IN+2 

GO TO 3324 

3321 CONTINUE 

DO 3322 K=1,NP 

3322 KM(IN+l,K)=KB(L+l,K)+KB(L+4,K)+KB(LL+l,K)+KB(LL+4 ,K) 
IN=IN+1 

3324 CONTINUE 

3325 CONTINUE 
IF(IJ.LT.l) GO TO 3345 

9085 F0RMAT(1H ,6F12.0) 

DO 3330 I“1,IJ 
DO 3330 J=l,7 
3330 FS(I,J)=0.0 

DO 3340 1=1,6 
DO 3335 J=1,IN 
DO 3335 K=1,NP 

3335 FS(J,7)=FS(J,7)+KM(J,K) *UVW(K,I) 

DO 3340 J=1,IJ 
FS(J,I)=FS(J,7) 

3340 FS(J,7)=0.0 
3345 CONTINUE 
C 

C USE DISPLACEMENT B.C.S TO ELIMINATE DISPLACEMENTS 
C 

INN=IN 

IN=0 

IF(IJ.LT.l) GO TO 3575 

IF( (NBX.LE. 1) .OR. (NBY.LE. 1) ) GO TO 3400 
DO 3375 1=2, NBX 
DO 3375 J=2,NBY 
L=( (I-l) * (NBY+1)+(J-1) ) *6 
DO 3350 K=1,INN 

KN(K, IN+l)=KM(K,L+l)+KM(K,L+4) 
KN(K,IN+2)=KM(K,L+2)+KM(K,L+5) 

3350 KN(K,IN+3)=KM(K,L+3)+KM(K,L+6) 

IN=IN+3 
3375 CONTINUE 
3400 CONTINUE 

IF(NBX.LE.l) GO TO 3475 
DO 3450 1=2, NBX 
L=6*(I-1)* (NBY+1) 

LL=NBY*6+L 
IF(LY) GO TO 3430 
DO 3425 K=1,INN 

KN(K,IN+l)=KM(K,L+l)+KM(K,L+4)+KM(K,LL+l)+KM(K,LL+4) 
KN (K , IN+2 ) =KM (K , L+2 ) +KM ( K , L+5 ) +KM (K, LL+2 ) +KM (K , LL+5 ) 
3425 KN(K,IN+3)=KM(K,L+3)+KM(K,L+6)+KM(K,LL+3)+KM(K,LL+6) 
IN=IN+3 
GO TO 3450 
3430 CONTINUE 
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IF(ISYM.EQ. 1) GO TO 3445 
DO 3440 K=1,INN 

KN(K, IN+1)=KM(K,L+1)+KM(K, L+4) 

KN(K, IN+2)=KM(K,L+3)+KM(K,L+6) 

KN ( K , IN+3 ) =*KM ( K , LL+1 ) +KM ( K , LL+4 ) 

3440 KN(K, IN+4)=KM(K,LL+3)+KM(K,LL+6) 

IN=IN+4 
GO TO 3449 
3445 CONTINUE 

DO 3447 K=l, INN 
KN(K,IN+l)=KM(K,L+2)+KM(K,L4-5) 

3447 KN(K,IN+2)=KM(K,LL+2)+KM(K,LL+5) 

IN=IN+2 

3449 CONTINUE 

3450 CONTINUE 
3475 CONTINUE 

IF(NBY.LE,1) GO TO 3510 
DO 3500 J=2,NBY 
L=6*(J-1) 

LL=6*NBX* (NBY+1) +L 
IF (LX) GO TO 3490 
DO 3480 K=1,INN 

KN(K,IN+l)=*KM(K,L+l)+KM(K,L+4)+KM(K,LL+l)+KM(K,LL+4) 
KN(K, IN+2)=KM(K,L+2)+KM(K,L+5)+KM(K,LL+2)+KM(K,LL+5) 
3480 KN(K,IN+3)=KM(K,L+3)+KM(K,L+6)+KM(K,LL+3)+KM(K,LL+6) 
IN=IN+3 
GO TO 3500 
3490 CONTINUE 

IF(ISYM.EQ. 1) GO TO 3496 

DO 3495 K=1,INN 

KN ( K , I N-f- 1 ) = KM ( K , L+ 2 ) + KM ( K , L+ 5 ) 

KN(K,IN+2)=KM(K,L+3)+KM(K,L+6) 

KN ( K , IN+3 ) ==KM ( K , LL+2 ) +KM ( K , LL+5 ) 

3495 KN(K,IN+4)=KM(K,LL+3)+KM(K,LL+6) 

IN=IN+4 

GO TO 3499 

3496 CONTINUE 

DO 3497 K=1,INN 
KN(K,IN+1)=KM(K,L+1)+KM(K, L+4) 

3497 KN(K,IN+2)=KM(K,LL+l)+KM(K,LL+4) 

IN=IN+2 

3499 CONTINUE 

3500 CONTINUE 
3510 CONTINUE 

IF(.NOT.LY) GO TO 3550 
L=6*NBY 

LL=L+6*NBX* (NBY+1) 

IF(ISYM.EQ. 1) GO TO 3530 
DO 3520 K=1,INN 

KN(K,IN+l)==KM(K,L-M)+KM(K,L+4)+KM(K,LL+l)+KM(K, LL+4) 
3520 KN(K,IN+2)=KM(K,L+3)+KM(K, L+6) +KM(K, LL+3 ) +KM (K, LL+6) 
IN=»IN+2 
GO TO 3540 
3530 CONTINUE 

DO 3535 K=l, INN 

3535 KN(K, IN+l)=KM(K,L+2)+KM(K,L+5)+KM(K,LL+2)+KM(K,LL+5) 
IN=IN+1 
3540 CONTINUE 
3550 CONTINUE 

IF(.NOT.LX) GO TO 3570 
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I>=6*NBX* (NBY+1) 

LL=L-^6*NBY 

IF(ISYM.EQ. 1) GO TO 3565 
DO 3560 K=1,INN 

KN ( K , I N + 1 ) = KM ( K , L+ 2 ) -f KM ( K , L+ 5 ) + KM ( K / LL+ 2 ) + KM ( K , LL+ 5 ) 
3 560 KN(K, IN+2)=KM(K,L-f3)+KM(K,L+6)+KM(K,LL+3)+KM(K,LL+6) 
IN=IN+2 
GO TO 3569 

3565 CONTINUE 

DO 3566 K=1,INN 

3566 KN{K,IN+l)=KM(K,L+l)+KM(K,L+4)+KM(K,LL+l)+KM(K,LL+4) 
IN=IN+1 

3569 CONTINUE 

3570 CONTINUE 
C 

C GET UNCONSTRAINED DISPLACEMENTS FOR 6 UNIT STRAIN CASES 
C 

12=0 

CALL GELIM(MMM,IJ,KN,7,FS,IP,IZ,WK,IERR) 

C 

C GET A COMPLETE SET OF TOTAL DISPLACEMENTS 
C 

3575 IN=0 

3590 IF( (NBX. LE. 1) -OR. (NBY.LE. 1) ) GO TO 3650 
DO 3625 1=2, NBX 
DO 3625 J=2,NBY 
LL=(I*NBY+I+J-NBY-2) *6 
DO 3600 K=l,3 
KP=K+3 

DO 3600 L=l,6 

UVW ( LL+K , L) =UVW ( LL+K , L) -FS ( IN+K , L) 

3600 UVW(LL+KP,L)=UVW(LL+KP,L) -FS (IN+K,L) 

IN=IN+3 
3625 CONTINUE 

3650 IF(NBX.LE.l) GO TO 3750 
DO 3740 1=2, NBX 
L=6* (I-l) * (NBY+1) 

LL=L+NBY*6 
IF(LY) GO TO 3710 
DO 3700 K=l,3 
KP=K+3 

DO 3700 M=l,6 

UVW(L+K,M) =UVW(L+K,M) -FS(IN+K,M) 

UW(L+KP,M) =UVW(L+KP,M) -FS(IN+K,M) 
UVW(LL+K,M)=UW(LL+K,M) -FS (IN+K,M) 

3700 UVW(LL+KP,M)=UVW(LL+KP,M)-FS(IN+K,M) 

IN=IN+3 
GO TO 3740 
3710 CONTINUE 

IF(ISYM.EQ. 1) GO TO 3730 
DO 3720 M=l,6 

UVW ( L+1 , M) =UVW ( L-4-1 , M) -FS ( IN+1 , M) 

UVW(L+3 ,M)=UVW(L+3,M) -FS (IN+2,M) 

UVW (L+4 , M) =UVW ( L+4 , M) -FS ( IN+1 , M) 
UVW(L+6,M)=UVW(L+6,M) -FS (IN+2 ,M) 

UVW ( LL+1 , M) =UVW ( LL+1 , M) -FS ( IN+3 , M) 

UVW (LL+3 , M) =UVW ( LL+3 , M) -FS ( IN+4 , M) 

UVW (LL+4 , M) =UVW (LL+4 , M) -FS ( IN+3 , M) 

3720 UVW(LL+6,M)=UVW(LL+6,M) -FS (IN+4,M) 

IN=IN+4 
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GO TO 3737 
3730 CONTINUE 

DO 3735 M=l,6 

UVW ( L+2 , M) =UVW ( L+2 , M) ~FS ( IN+1 , M) 

UVW ( L+5 , M) =UVW ( L+5 , M) -FS ( IN+1 , M) 

UVW ( LL+2 M) =UVW ( LL+2 , M) -FS ( IN+2 , M) 

373 5 UW(LL+5,M) =UVW(LL+5,M) ~FS(IN+2 ,M) 
IN=IN+2 
3737 CONTINUE 
3740 CONTINUE 

3750 IF(NBY.LE.l) GO TO 3850 
DO 3845 J=2,NBY 
D= 6*(J-1) 

LL=L-»- 6* (NBY+1) *NBX 
IF(LX) GO TO 3810 
DO 3800 K=l,3 
KP=K+3 

DO 3800 M“l,6 

UVW(L+K,M) =UVW(L+K,M) -FS(IN-fK,M) 
UVW(L+KP,M) =UVW(L+KP,M) -FS(IN+K,M) 
UVW(LL+K,M) =UVW(LL+K,M) -FS(IN+K,M) 
3800 UVW(LL+KP,M) =UVW (LL+KP,M) -FS(IN+K,M) 
IN=«IN+3 
GO TO 3845 
3810 CONTINUE 

IF{ISYM.EQ. 1) GO TO 3825 
DO 3820 M=l, 6 

UVW ( L4-2 , M ) ==UVW ( L+2 , M) -FS ( IN+ 1 , M ) 
UVW(L+3,M)=UVW(L+3,M) -FS(IN+2,M) 
UVW(L+5,M)=UVW(L+5,M) -FS (IN+1,M) 
UVW(L+6,M)=UVW(L+6,M) -FS(IN+2,M) 

UVW ( LL+2 , M) ==UVW (LL+2 , M) -FS ( IN+3 , M) 
UVW(LL+3 ,M)=UVW(LL+3 ,M) -FS(IN+4,M) 
UVW(LL+5,M)=UVWCLL+5,M) -FS (IN+3,M) 

3820 UVW(LL+6,M)=UVW(LL+6,M) -FS (IN+4 ,M) 
IN=IN+4 
GO TO 3835 
3825 CONTINUE 

DO 3830 M=l,6 

UVW ( L+ 1 , M) =UVW ( L+l , M) -FS ( IN+ 1 , M) 
UVW(L+4 ,M)=UVW(L+4 ,M) -FS(IN+1,M) 
UVW(LL+1,M)=UVW(LL+1,M) -FS(IN+2,M) 

3830 UVW(LL+4,M)=UVW(LL+4,M)-FS(IN+2,M) 
IN=IN+2 
3835 CONTINUE 
3845 CONTINUE 
3850 CONTINUE 

IF(.NOT.LY) GO TO 3870 
L=6*NBY 

LL=E+6*NBX* (NBY+1) 

IF(ISYM.EQ.l) GO TO 3865 
DO 3860 M=l,6 

UVW(L+1,M)=UVW(L+1,M) -FS (IN+1,M) 

UVW ( L+3 , M) =UVW ( L+3 , M) -FS ( IN+2 , M) 

UVW ( L+4 , M) -UVW ( L+4 , M) -FS ( IN+ 1 , M) 
UVW(L+6,M)=UVW(L+6,M) -FS (IN+2,M) 

UVW ( LL+1 , M) =UVW ( LIi+1 , M) -FS ( IN+1 , M) 
UVW(LL+3,M)=UVW(LL+3,M) -FS(IN+2,M) 

UVW (LL+4 , M) =UVW (LL+4 , M) -FS ( IN+1 , M) 

3860 UVW(LL+6,M)=UVW(LL+6,M)-FS(IN+2,M) 
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non o n o 


IN=IN+2 
GO TO 3869 

3865 CONTINUE 

DO 3866 M=l,6 

UVW ( L+2 , M) =UVW ( L+2 , M) -FS ( IN+1 , M) 

UVW (L+5 , M) =»UVW ( L+5 , M) -FS ( IN+1 , M) 

UVW ( LL+2 , M) -UVW ( LL+2 , M) -FS ( IN+1 , M) 

3866 UVW(LL+5,M)=UVW(LL+5,M) -FS (IN+1,M) 

IN=IN+1 

3869 CONTINUE 

3870 CONTINUE 

IF (.NOT. LX) GO TO 3890 
L=6*NBX*(NBY+1) 

LL=L+6*NBY 

IF(ISYM.EQ. 1) GO TO 3885 
DO 3880 M“l,6 

UVW ( L+2 , M) =UVW ( L+2 , M) -FS ( IN+1 , M) 

UVW ( L+3 , M) =UVW ( L+3 , M) *-FS ( IN+2 , M) 
UVW(L+5,M)=UVW(L+5,M)-FS(IN+1,M) 
UVW(L+6,M)=UVW(L+6,M) -FS (IN+2 ,M) 

UVW { LL+2 , M) =UW ( LL+2 , M) -FS ( lN+1 , M) 

UVW ( LL+3 , M) =UVW ( LL+3 , M) -FS ( IN+2 , M) 
UVW(LL+5,M)=UVW(LL+5,M) -FS (IN+1,M) 

3880 UVW(LL+6,M) =UVW(LL+6,M) -FS (IN+2 ,M) 

IN=IN+2 
GO TO 3889 

3885 CONTINUE 

DO 3886 M=l,6 

UVW(L+1,M)=UW(L+1,M) -FS (IN+1,M) 

UVW(L+4,M)=UVW(L+4 ,M) -FS (IN+1,M) 

UVW(LL+1,M) =UVW(LL+1,M) -FS (IN+1,M) 

3886 UVW(LL+4,M)=UVW(LL+4,M) -FS(IN+1,M) 

IN=IN+1 

3889 CONTINUE 

3890 CONTINUE 

COMPUTE NODAL FORCES 

IF(NP.LT.l) GO TO 3950 
DO 3900 1=1, NP 
DO 3900 J=l,6 
FB(I,J)=0.0 
DO 3900 K=1,NP 

3900 FB(I,J)=FB(I, J)+KB(I,K) *UVW(K, J) 

3950 CONTINUE 

COMPUTE SIDE LOADS FOR EACH OF THE 6 UNIT STRAIN CASES 

YLC=YL*C 
XLC=XL*C 
DO 4500 1=1,6 
DO 4200 J=1,NBY+1 
K=(J-1)*6 

FT(l,I)=FT(l,I)-4-FB(K+l,I)+FB(K+4,I) 
FT(4,I)=FT(4,I)+FB(K+2, I)+FB(K+5,I) 

4200 FT(6,I)=FT(6,I)+FB(K+3,I)+FB(K+6,I) 

FT ( 1 , I ) =-FT ( 1 , I ) / YLC 
FT(4,I)=-FT(4,I)/YLC 
FT ( 6 , I ) =-FT ( 6 , I ) / YLC 
DO 4300 J=1,NBX+1 
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o o o 


4300 


K=(NBY+1) * (J-1) *6 

FT(2,I)==FT(2,I)+FB(K+2,I)+FB(K+5,I) 
FT(5,I)«FT(5,I)-fFB(K+3,I)+FB(K+6,I) 
FT(2,I)=-FT(2,I)/XLC 
FT ( 5 , I ) =-FT ( 5 , I ) /XLC 
DO 4 4 00 J=1,NBX4-1 
DO 4400 K=1,NBY+1 
L=( (NBY+1) * (J-1)+(K“1) ) *6 
4400 FT(3,I)=FT(3 ,I)+FB(L+3,I) 

FT(3,I)=“FT{3,I)/ (XL*YL) 

4500 CONTINUE 

IF(ISYM.EQ. 0) THEN 
DO 4600 1=1,6 
DO 4600 J=l,6 
TF(I, J)=FT(I,J) 

4600 FT(I,J)=0.0 
END IF 

IF( (LX.OR.LY) .AND. (ISYM.EQ.O) ) GO TO 2440 
IF(ISYM.EQ. 1) THEN 
DO 4700 1=1,6 
TF(I,4)=FT(I,4) 

IF(LY) TF(I,5)=FT(I,5) 

4700 IF(LX) TF(I,6)=FT(I, 6) 

END IF 

DO 4800 1=1, 6 
DO 4800 J=l,6 
4800 FT(I, J)=TF(I, J) 

CALCULATE THE ELASTIC CONSTANTS OF THE UNIT CELL 


CALL MATO PS ( KA , FT , DET , HH ) 
EX=1.0/FT(1, 1) 
EY=1.0/FT(2,2) 

EZ=1. 0/FT(3 ,3) 

GXY=1. 0/FT(4 , 4 ) 
GY2=1.0/FT(5, 5) 
GXZ=1.0/FT(6, 6) 

PRXY=-FT (1,2)/FT(2,2) 
PRXZ=-FT (1,3) /FT (3,3) 
PRYZ=-FT (2,3) /FT (3,3) 
MIXXY=FT(4, 1)/FT(1, 1) 
MIXXZ=FT(6, 1)/FT(1,1) 
MIXYZ=FT(5, 1)/FT(1, 1) 
MIYXY=FT(4, 2)/FT(2,2) 

MI YXZ=FT (6,2) /FT (2,2) 
MIYYZ=FT (5,2) /FT (2,2) 
MIZXY=FT(4, 3)/FT(3, 3) 
MIZXZ=FT(6, 3)/FT(3, 3) 
MIZYZ=FT (5,3) /FT (3,3) 
CXYXZ=FT(4 ,5)/FT(5,5) 
CXYYZ=FT(4, 6)/FT(6,6) 
CXZYZ=FT(5,6)/FT(6, 6) 
WRITE (6, 9560) 

WRITE (6, 9560) 

WRITE (6, 9600) 

WRITE (6, 9560) 

WRITE (6, 9500) 

WRITE (6, 95 10) 

WRITE (6, 9520) 

WRITE (6, 9530) 


EX,EY,EZ 
GYZ , GXZ , GXY 
PRYZ , PRXZ , PRXY 
MIXYZ,MIYY2,MIZYZ 
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WRITE (6, 9540) 
WRITE(6,9550) 
WRITE(6,9550) 
WRITE (6, 9 560) 
WRITE (6, 9560) 
9000 FORMAT(F12.5) 
9005 FORMAT (IH , 8F9 
9010 FORMAT(F12. 2) 
9020 FORMAT(6F12 .2) 

9030 FORMAT(I5) 

9031 FORMAT(3I5) 
9040 FORMAT (L5) 

9050 FORMAT (IH , 
9055 FORMAT (IH , 
9060 FORMAT (IH , 
9070 FORMAT (IH , 
9080 FORMAT (IH , 
9090 FORMAT (IH , 
9100 FORMAT (IH , 
9110 FORMAT (IH , 
9120 F0RMAT(1H , 
9130 FORMAT(lH , 
9140 FORMAT (IH , 
9150 FORMAT (IH , 
9160 FORMAT (IH , 
9170 FORMAT (IH , 
9180 FORMAT (IH , 
9190 FORMAT (IH , 
9200 FORMAT (IH , 
9210 FORMAT (IH , 
9220 FORMAT(lH , 
9230 FORMAT (IH , 
9240 FORMAT (IH , 
9250 F0RMAT(1H , 
9260 FORMAT (IH , 
9270 FORMAT (IH , 
9280 FORMAT (IH , 
9290 FORMAT (IH , 
9300 FORMAT (IH , 
9310 FORMAT (IH , 
9320 FORMAT (IH , 
9330 FORMAT (IH , 
9390 F0RMAT(1H , 
9400 FORMAT (IH , 
9420 FORMAT (IH , 

9430 F0RMAT(1H , 
9440 F0RMAT(1H , 

9450 FORMAT (IH , 

9460 FORMAT (IH , 

9470 FORMAT (IH , 

9480 FORMAT (IH , 

9490 FORMAT (IH , 

9500 FORMAT (IH , 

9510 FORMAT (IH , 

9520 FORMATdH , 

9530 FORMAT (IH , 

9540 FORMAT (IH , 

9550 FORMAT (IH , 

9560 FORMAT (IH ) 

9600 FORMAT (IH 


MIXXZ , MIYXZ , MIZX2 
MIXXY,MIYXY,MIZXY 
MIXXY , MIYXY , MIZXY 


3) 


SYMMETRY ABOUT Y2 PLANE (T OR F) ') 

SYMMETRY ABOUT XZ PLANE (T OR F) ' ) 

INPUT SUBCELL ID. NO, AT LOCATION •, 214 ) 
INPUT SUBCELL ROTATION CODE NO. AT ',2 14) 
INPUT NO. SUBCELLS (X DlR. ) IN UNIT CELL') 
INPUT NO. SUBCELLS (Y DIR.) IN UNIT CELL’) 
INPUT NO. COMPOSITE MATERIALS NEEDED, NM') 
INPUT NO. OF MASTER SUBCELLS NEEDED, NB') 
INPUT E IN FIBER DIRECTION*) 

INPUT E NORMAL TO FIBER DIRECTION') 

INPUT MAJOR POISSONS RATIO IN LT PLANE ' ) 
INPUT POISSONS RATIO IN TT PLANE') 

INPUT SHEAR MODULUS G IN LT PLANE ’ ) 

INPUT SHEAR MODULUS G IN TT PLANE') 

SELECT A MATERIAL NUMBER FROM ONE TO TEN') 
MATERIAL PROPERTY DATA ECHO ' ) 

GRID GEOMETRY FOR SUBCELL NUMBER ’,15) 
INPUT SIDE LENGTH OF SUBCELL (X DIR.) •) 
INPUT SIDE LENGTH OF SUBCELL (Y DIR.) ') 
INPUT HEIGHT OF SUBCELL (Z DIR.) ') 

INPUT NO. INTEG/PTS. ON SUBCELL (X DIR.)’) 
INPUT NO. INTEG/PTS. ON SUBCELL (Y DIR.)') 
INPUT NO. INTEG/PTS, ON SUBCELL (Z DIR,)') 
INPUT DELTA X IN % FOR INTERVAL NO. ’,15) 
INPUT DELTA Y IN % FOR INTERVAL NO. *,I5) 

FOR INTERVAL NO. ’,15) 


Z IN % 


,314) 


INPUT DELTA 

INPUT NO. OF MATLS.IN SUBCELL NO. ’,15) 
SPECIFY THE CURRENT MAIL. ID. NO.') 

SUBCELL HT. NOT EQUAL UNIT CELL HT. ') 

INPUT 1ST SPH. FIBER ANG. AT GRID PT. 

INPUT 2ND SPH. FIBER ANG. AT GRID PT.',3I4) 
INPUT OCTANT NO. AT GRID NO. ',314) 

INPUT MATL. VOL/FR. AT OCTANT NO. ',15) 

INPUT SIDE LENGTH OF UNIT CELL IN X DIR.') 
INPUT SIDE LENGTH OF UNIT CELL IN Y DIR.') 
INPUT DIST.(%) ORIGIN TO UNIT CELL NODE ' , 13 ) 
INPUT THICKNESS OF UNIT CELL IN INCHES') 
INPUT 1ST FIBER SPHERICAL ANGLE') 

INPUT 2ND FIBER SPHERICAL ANGLE ' ) 

EX,EY,EZ = ' ,3F12 .2) 

GYZ,GXZ,GXY = ',3F12.2) 


MUZY,MUZY,MUYX = 
NUYZ,X ; NUYZ,Y 
NUXZ , X ; NUXZ , Y 
NUXY,X ; NUXY,Y 


, 3F12 .4) 
NUYZ,Z = 
NUXZ,Z = 
NUXY,Z = 


, 3F12 
, 3F12 


4) 

4) 


' ,3F12.4) 


, 13 X, 'ELASTIC CONSTANTS OF THE COMPOSITE ') 
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non 


5000 END 

SUBROUTINE TO TRANSFORM THE SUBCELL STIFFNESS MATRIX 

SUBROUTINE TRNS(QT,KS) 

DIMENSION SK(24,24) 

REAL KS (24 , 24) 

INTEGER QT 

7117 FORMAT (IH , ' IN SUB TRNS , QT=',I3) 

IFCQT.GE. 100) THEN 
DO 10 1=1,3 
SN=1. 0 

IF (I.EQ.3) SN=-1.0 
DO 10 J=l,24 
SK(I,J)=KS(I+3, J)*SN 
SK(I+3 , J) =KS (I , J) *SN 
SK(I+6, J)=KS (1+9, J) *SN 
SK(I+9, J)=KS (1+6, J) *SN 
SK (1+12, J)=KS (1+15, J) *SN 
SK(I+15,J)=KS(I+12,J) *SN 
SK(I+18, J)=KS (1+21, J) *SN 
10 SK(I+21, J)=KS (1+18, J) *SN 
DO 20 J=l,3 
SN=1. 0 

IF(J.EQ.3) SN=-1.0 
DO 20 1=1,24 
KS(I,J)=SK(I, J+3) *SN 
KS (I, J+3)=SK(I, J) *SN 
KS(I, J+6)=SK(I,J+9) *SN 
KS(I, J+9)=SK(I, J+6) *SN 
KS(I, J+12)=SK(I,J+15) *SN 
KS(I, J+15)=SK(I,J+12) *SN 
KS(I,J+18)=SK(I,J+21) *SN 
20 KS (I, J+21)»SK(I, J+18) *SN 
QT=QT“100 
END IF 

IF(QT.GE.IO) THEN 
DO 100 1=1,3 
SN=-1.0 

IF(I.EQ. 1) SN=1.0 
DO 100 J=l,24 
SK(I, J)=KS(I+9,J) *SN 
SK(I+3,J)=KS(I+6,J)*SN 
SK(I+6, J)=KS (1+3 , J) *SN 
SK(I+9, J)=KS(I, J) *SN 
SK(I+12, J)=KS (1+21, J) *SN 
SK(I+15, J)=KS(I+18, J) *SN 
SK ( 1+18, J)=KS( 1+15, J) *SN 
100 SK(I+21, J)=KS(I+12, J) *SN 
DO 200 J=l,3 
SN=-1.0 

IF(J.EQ,1) SN=1.0 
DO 200 1=1,24 
KS(I,J)=SK(I, J+9) *SN 
KS(I, J+3)=SK(I,J+6) *SN 
KS(I, J+6)=SK(I,J+3) *SN 
KS(I, J+9)=SK(I, J) *SN 
KS(I, J+12)=SK(I,J+21) *SN 
KS(I, J+15)=SK(I,J+18) *SN 
KS(I, J+18)=SK(I,J+15) *SN 
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200 KS (I, J+21)=SK(I, J+12) *SN 
QT=QT-10 
END IF 

IF( (QT.EQ.l) .OR. (QT.EQ. 3) ) THEN 

DO 300 1=1,3 

SN=1.0 

IF(I.EQ.2) SN=~1.0 
IF(I.EQ.l) 11=1+1 
IF(I.EQ.2) 11=1-1 
IF(I.EQ. 3) II=I 
DO 300 J=l, 24 
SK(I , J) =KS (11+12 ,J) *SN 
SK(I+3 , J)=KS(II+15, J) *SN 
SK(I+6 ,J)=KS(II ,J)*SN 
SK(I+9 ,J)=KS(II+3 ,J)*SN 
SK ( 1+12 ,J)=KS (11+18, J) *SN 
SK(I+15, J)=KS (11+21, J) *SN 
SK(I+18 , J)=KS (II+6 ,J)*SN 
300 SK(I+21, J)=KS(II+9 , J) *SN 
DO 400 J=l,3 
SN=1.0 

IF(J.EQ.2) SN=-1.0 
IF(J.EQ.l) JJ=J+1 
IF(J.EQ.2) JJ=J-1 
IF(J.EQ.3) JJ=J 
DO 400 1=1,24 
KS(I,J )=SK(I, JJ+12) *SN 
KS(I,J+3 )=SK(I, JJ+15) *SN 
KS(I,J+6 )=SK(I,JJ )*SN 
KS(I,J+9 )=SK(I,JJ+3 )*SN 
KS(I,J+12)=SK(I, JJ+18) *SN 
KS(I, J+15)=SK(I, JJ+21) *SN 
KS(I,J+18)=SK(I,JJ+6 )*SN 
400 KS (I, J+21)=SK(I, JJ+9 ) *SN 
END IF 

IF( (QT.EQ. 2) .OR. (QT.EQ. 3) ) THEN 
DO 500 1=1,3 
SN=-1. 0 

IF(I.EQ.3) SN=1.0 
DO 500 J=l,24 
SK(I, J)=KS (1+18 , J) *SN 
SK(I+3 , J)=KS (1+21, J) *SN 
SK(I+6, J)=KS(I+12,J) *SN 
SK(I^9 , J)=KS (1+15, J) *SN 
SK ( 1+12 ,J) =KS (1+6, J) *SN 
SK(I+15, J)=KS(I+9,J) *SN 
SK(I^18, J)=KS(I, J) *SN 
500 SK(I-^21,J)=KS(I + 3,J)*SN 
DO 600 J=l,3 
SN=-1. 0 

IF(J.EQ.3) SN=1.0 
DO 600 1=1,24 
KS(I,J)=SK(I,J+18) *SN 
KS (I, J+3)=SK(I, J+21) *SN 
KS(I, J+6)=SK(I, J+12) *SN 
KS (I, J+9)»SK(I, J+15) *SN 
KS (I, J+12)=SK(I, J+6) *SN 
KS(I,J+15)=SK(I,J+9)*SN 
KS (I,J+18)=SK(I, J) *SN 
600 KS(I, J+21)=SK(I, J+3) *SN 
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END IF 
900 RETURN 
END 
C 

C MATRIX INVERSION 
C 

SUBROUTINE INV(AK) 

DIMENSION KA(7) ,AK(9,9) 

KA(1)=10 

KA(2)=9 

KA(3)=9 

KA(4)>=1 

KA(5)=9 

KA(6) =0 

KA(7)=0 

CALL MATOPS (KA,AK,DET, DUMMY) 

RETURN 

END 

/ 
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A££EMDIK-IL.,- ...I NP UT DAT A . F ORM AT 

This program operates interactively and is largely self-explanatory. The first 
block of input data establishes a 2-D array containing principal elastic constants 
for each different material needed in the run. Each material is assumed to be 
orthotropic with a principal plane of isotropy. Therefore, only six elastic 
constants {per material) must be specified. The first input is an integer between 
one and ten quantifying the number of materials needed for the run. This number 
sizes the array that contains the elastic constants of each material. For 
convenience six sets of elastic constants reside in the program in the form of 
data statements: 

1. High Modulus Graphite/Epoxy (E^ = 30 msi) 

2. Intermediate Modulus Graphite/Epoxy (E^ *» 25 msi) 

3. Low Modulus Graphite/Epoxy (E^ = 20 msi) 

4. Fiberglass/Epoxy (E^ = 10 msi) 

5. Aluminum {E = 10 msi) 

6. Bulk Epoxy (E = . 5 msi) 

When called for, by integer 1-6, the corresponding set of elastic constants is 
inserted into the next empty row of the material property array. The material 
calling sequence establishes a new material numbering sequence for subsequent use. 
Any material may be requested more than once. If a material number between 7 and 
10, inclusive, is requested the program immediately asks for six elastic constants 
to be input. Young's modulus in the reinforcement direction is requested first, 
followed by Young's modulus in the plane of isotropy. Then the program requests 
two Poisson's ratios: first, that which specifies contraction normal to fiber 

direction (per unit strain in the fiber direction); second, that which relates to 
the plane of isotropy. Finally, two shear moduli are requested: first, in a plane 
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containing the fiber principal axis; second, in the plane of isotropy. The 
material sequencing then resumes. The number of materials listed must equal the 
initial input integer specifying the number of materials to be used. 

The next number requested (an integer from one to ten) establishes the total 
number of master subcells needed to build the unit cell. 

The next block of input data establishes size, subcell division, and subcell 
content of the unit cell model. This data describes the unit cell as a 2-D 
compartmented container for subcells. The first number (an integer from one to 
ten) specifies the number of subcell compartments in the unit cell, as the cell is 
traversed in the x-direction (Figure 38) . The second integer specifies the number 
of subcell compartments in the y direction. The program assumes single subcell 
thickness in the z-direction of the unit cell. (Other versions of this program 
allow more than one subcell in the thickness direction.) 

The next series of real numbers establish dimensions of all unit cell 
compartments. The first number sets unit cell height, i.e., ply thickness and 
subcell height. The second real number sets the x-parallel side length of the 
unit cell. This number is followd by a sequence of percentages (0.0 to 100.0) 
specifying the fractions of x-parallel side length from the origin to the largest 
x-dimension within each subcell (Figure 38) . 

The next set of real numbers establishes the same dimensions along a y-parallel 
side of the unit cell. Numbers and sizes of all unit cell compartments are now 
fixed. The next pair of integers specifies both the type of master subcell 
associated with the first compartment of the unit cell and the type of coordinate 
transformation to be performed on that master subcell (prior to placement in that 
unit cell compartment) . The first integer (specifying subcell type) also tacitly 
establishes the sequence of master subcell descriptions that follow. The second 
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integer selects the desired subcell coordinate transformation according to the 
following code: 


Integer 

Transformation (When viewed from the subcell origin, 
looking in the + direction along a coordinate axis) 

-1 

No subcell inserted 

0 

No subcell transformation necessary 

1 

Rotation Of 90^ CCW about z axis 

2 

Rotation of 180° CCW about z axis 

3 

Rotation of 270® CCW about z axis 

IX 

Rotation of 270° CCW about x axis (plus whichever one 


of the previous transformations is specified by X) 

IXX 

Reversal of the z axis (plus whichever one of the 


previous transformations is specified by XX) 


CW Indicates clockwise transformation 

CCW Indicates counterclockwise transformation 


Subsequent pairs of integers are coded in the order shown in Figure 39 until all 
the compartments of the unit cell are filled with transformed subcells. There is 
no way to put more than one subcell in a compartment. This completes the 
description of the unit cell. 

It remains to describe the details of material placement within each master 
subcell. To facilitate this task an outer do-loop is established over the number 
of master subcells. Each master subcell stiffness matrix is then formed by an 
inner do-loop ranging over the number of materials in that subcell. The volume 
integration associated with JJJ [B]"^ [D] [B] d(vol) is performed separately for each 
material. The same integration scheme and 3-D integration network is used for 
each material in the subcell. However, different integration networks can be used 
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for different master subcells. What follows pertains to the integration grid 
dimensions for the first master subcell. 

The next real number input is the integration grids dimension in the x-direction, 
followed by an integer that specifies the number of integration points on the x- 
axis. This is followed by a sequence of percentages (0.0 to 100.00) specifying 
the fraction of the overall grid dimension represented by the distance of each x- 
axis integration point from the origin (except for the first and last values which 
are assumed to be 0.0% and 100.0%). 

The next two batches of input data repeat this same information for the 
integration grid points along the y-axis and along the z-axis. Subsequently, an 
integer is specified establishing the number of materials in the first master 
subcell. This sets the upper limit on the inner do-loop that begins at this 
point. 

The first integer related to the inner do-loop is a material number. This selects 

a set of material properties from that row of the material property matrix. The 

program is now ready for detailed information on fiber orientation and amount of 
material at each integration point. First, two spherical angles 

degrees) specifying fiber direction (Figure 11) at the first integration point, 
are input, followed by the volume fraction of each octant of volume surrounding 
that integration point that contains the material of interest. These volume 
fractions are determined by examining photomicrographs of the actual composite and 
unit cell models constructed from the same photomicrographs. To save time each 
octant volume fraction is read in by means of a two number code. The first number 
specifies which of the eight octants is being referenced. Octant one is the +++ 
octant, octant two the ++- one, octant three the +-+ one, octant four the + — one, 
octant five the -++ one, octant six the -+- one, octant seven the one, and 
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octant eight the 


one (Figure 12) . A negative octant reference means skip to 

the next integration point. Any octant reference number between one and eight 
must be followed by a real number (between 0.0 and 1.0) specifying the volume 
fraction of current material occupying that octant of that integration point. 

The sequence of points over which the intregration ranges is fixed. The first 
point is always the origin. The path then runs along the plus z-axis, returns to 
the next integration point along the y-axis, and moves in the plus z-direction 
again, etc. (Figure 40) . When the farthest integration point from the origin is 
reached the program returns to the start of the inner do loop for the inclusion of 
the next material within that master subcell. When all materials in the master 
subcell have been allotted to the various integration points, the program returns 
to the outer do loop to repeat the same input data pattern for the next master 
subcell. When all the master subcells have been described the program proceeds to 
calculate the unit cell moduli. 
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APPENDIX III 


Sample Problem 

As an example of data input/output consider the previous model of the plain weave 
microgeometry based on diamond shaped tow cross-sections (Figure 24) . Only one 
master subcell is required. The unit cell consists of 16 inhomogeneous finite 
elements. Their eight-node stiffness matrices are all transformations of the 
master subcell stiffness matrix. For simplicity the master subcell may be 
considered to be a unit cube containing 50% bulk matrix material and 50% 
unidirectional composite tow material. 

The unit cell, before the subcells are inserted, may be considered to be an empty, 
four by four, compartmented container (Figure 24) , The first subcell is placed in 
the first compartment bounded by the planes x=0, x=l, y=0, y=l . The first 
subcells is the untransformed master subcell. The second subcell is placed in the 
second compartment bounded by the planes x=0, x«l, y=l, y=2 . It is a transformed 
master subcell. The transformation consists of a rotation of 90° about an axis 
normal to the middle plane of the fabric plus a reflections about the middle 
plane. The third subcell is placed in compartment x=0, x=l, y=2, y=3. It is a 
master subcell reflected about the fabric middle plane; etc. 

The following sequence of prompts and inputs leads to the composite moduli 
predictions that follow. All input values are enclosed in double brackets. Most 
of the master subcell inputs come from Table 2. 
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INPUT NO. COMPOSITE MATERIALS NEEDED, NM 
((3)) 

SELECT A MATERIAL NUMBER FROM ONE TO TEN 
((3)) 

SELECT A MATERIAL NUMBER FROM ONE TO TEN 
((3)) 

SELECT A MATERIAL NUMBER FROM ONE TO TEN 

(( 6 )) 

MATERIAL PROPERTY DATA ECHO 


20000000.00 

1500000.00 

.30 

.45 

700000.00 

700000.00 

20000000.00 

1500000.00 

.30 

.45 

700000.00 

700000.00 

500000.00 

500000.00 

.35 

.35 

185000.00 

185000.00 


INPUT NO. OF MASTER SUBCELLS NEEDED, NB 

{(!)) 

INPUT NO. SUBCELLS (X DIR.) IN UNIT CELL 
((4)) 

INPUT NO. SUBCELLS (Y DIR.) IN UNIT CELL 
((4)) 

INPUT THICKNESS OF UNIT CELL IN INCHES 

<( 1 . 0 )) 

INPUT THICKNESS OF UNIT CELL IN X DIRECTION 
( (4.0) ) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 2 
((25.0)) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 3 
( (50.0) ) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 4 
<<75. 0)) 

INPUT SIDE LENGTH OF UNIT CELL IN Y DIRECTION 
((4.0)) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 2 
((25.0)) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 3 
((50.0)) 

INPUT DIST. (%) ORIGIN TO UNIT CELL NODE 4 
((75.0)) 

INPUT SUBCELL ID. NO. AT LOCATION 1 1 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 1 1 

(( 0 )) 

INPUT SUBCELL ID. NO. AT LOCATION 1 2 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT I 2 
(dOD) 

INPUT SUBCELL ID. NO. AT LOCATION 1 3 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 1 3 

(( 100 )) 

INPUT SUBCELL ID. NO. AT LOCATION 1 4 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 1 4 

((D) 

INPUT SUBCELL ID. NO. AT LOCATION 2 1 

((D) 
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INPUT SUBCELL ROTATION CODE NO. AT 2 
( (103) ) 

INPUT SUBCELL ID. NO. AT LOCATION 2 

((!)) 

INPUT SUBCELL ROTATION CODE NO. AT 2 

(( 2 )) 

INPUT SUBCELL ID. NO. AT LOCATION 2 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 2 
((3)) 

INPUT SUBCELL ID. NO. AT LOCATION 2 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 2 

(( 102 )) 

INPUT SUBCELL ID. NO. AT LOCATION 3 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 3 

(( 100 )) 

INPUT SUBCELL ID. NO. AT LOCATION 3 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 3 

((D) 

INPUT SUBCELL ID. NO. AT LOCATION 3 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 3 

(( 0 )) 

INPUT SUBCELL ID. NO. AT LOCATION 3 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 3 
( (lOD) 

INPUT SUBCELL ID. NO. AT LOCATION 4 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 4 
((3)) 

INPUT SUBCELL ID. NO. AT LOCATION 4 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 4 

(( 102 )) 

INPUT SUBCELL ID. NO. AT LOCATION 4 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 4 
((103)) 

INPUT SUBCELL ID. NO. AT LOCATION 4 

((D) 

INPUT SUBCELL ROTATION CODE NO. AT 4 

(( 2 )) 

SYMMETRY ABOUT YZ PLANE (T OR F) 

( (F) ) 

SYMMETRY ABOUT XZ PLANE (T OR F) 

((F)) 


1 

2 

2 

3 

3 

4 
4 
1 
1 
2 
2 
3 

3 

4 
4 
1 
1 
2 
2 
3 

3 

4 
4 
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« ★ * 


GRID GEOMETRY FOR SUBCELL NUMBER 1 

* * ★ 

INPUT SIDE LENGTH OF SUBCELL (X DIR.) 

(( 1 . 0 )) 

INPUT NO. INTEG/PTS. OF SUBCELL (X DIR.) 

(( 2 )) 

INPUT SIDE LENGTH OF SUBCELL (Y DIR) 

(d.O)) 

INPUT NO INTEG/PTS. ON SUBCELL (Y DIR.) 

(( 2 )) 

INPUT HEIGHT OF SUBCELL (Z DIR.) 

(( 1 . 0 )) 

INPUT NO INTEG/PTS. ON SUBCELL (Z DIR.) 

(( 2 )) 

INPUT NO. OF MATLS. IN SUBCELL NO. 1 
((3)) 

SPECIFY THE CURRENT MATL. ID. NO. 

((D) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. Ill 

((D) 

INPUT MATL. VOL/FR AT OCTANT NO. 1 1 1 

((0.04)) 

INPUT OCTANT NO. AT GRID NO. Ill 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 112 

(( 2 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 1 12 

( ( 0 . 21 ) ) 

INPUT OCTANT NO. AT GRID NO. 112 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 1 2 1 

((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
( (90.0) ) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(dO.O)) 

INPUT OCTANT NO. AT GRID NO. 122 
((4)) 

INPUT MATL. VOL/FR AT OCTANT NO. 4 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 122 
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((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 211 
((5)) 

INPUT MATL. VOL/FR AT OCTANT NO. 5 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 211 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 212 

(( 6 )) 

INPUT MATL. yOL/FR AT OCTANT NO. 6 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 212 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 221 

((7)) 

INPUT MATL. VOL/FR AT OCTANT NO. 7 
( (0.04) ) 

INPUT OCTANT NO. AT GRID NO. 221 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 2 , 2 2 

(( 8 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 8 
((0.71)) 

INPUT OCTANT NO. AT GRID NO. 222 
((-!)) 

SPECIFY THE CURRENT MATL ID. NO. 

(( 2 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. Ill 

((D) 

INPUT MATL. VOL/FR AT OCTANT NO. 1 

(( 0 . 21 )) 

INPUT OCTANT NO. AT GRID NO. Ill 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
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(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 112 

(( 2 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 2 
((0.04)) 

INPUT OCTANT NO. AT GRID NO. 112 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 121 
((3) 

INPUT MATL. VOL/FR AT OCTANT NO. 3 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 121 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 122 
((4)) 

INPUT MATL. VOL/FR AT OCTANT NO. 4 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 122 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 211 
((5)) 

INPUT MATL. VOL/FR AT OCTANT NO. 5 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 21 1 

((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 212 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 221 
((7) 

INPUT MATL. VOL/FR AT OCTANT NO. 7 
((0.71)) 

INPUT OCTANT NO. AT GRID NO. 221 
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((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

((- 10 . 0 )) 

INPUT OCTANT NO, AT GRID NO. 22 2 

(( 8 )) 

INPUT MATL, VOL/FR AT OCTANT NO. 8 
((0,04)) 

INPUT OCTANT NO. AT GRID NO. 2 2 2 

((-!)) 

SPECIFY THE CURRENT MATL. ID. NO. 

((3)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. Ill 

((D) 

INPUT MATL. VOL/FR AT OCTANT NO. 1 
((0.75)) 

INPUT OCTANT NO. AT GRID NO. Ill 

((-D) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 1 12 

(( 2 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 2 
((0.75)) 

INPUT OCTANT NO. AT GRID NO. 112 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 121 
((3)) 

INPUT MATL. VOL/FR AT OCTANT NO. 3 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 12 1 

((-D) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 122 
((4)) 

INPUT MATL. VOL/FR AT OCTANT NO. 4 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 122 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 
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<( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO, 211 
((5)) 

INPUT MATL. VOL/FR AT OCTANT NO. 5 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 211 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 212 

(( 6 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 6 
((0.50)) 

INPUT OCTANT NO. AT GRID NO. 212 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 221 
((7)) 

INPUT MATL. VOL/FR AT OCTANT NO. 7 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 221 
((-!)) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT OCTANT NO. AT GRID NO. 222 

(( 8 )) 

INPUT MATL. VOL/FR AT OCTANT NO. 8 
((0.25)) 

INPUT OCTANT NO. AT GRID NO. 22 2 

((-!)) 


ELASTIC CONSTANTS OF THE COMPOSITE 


EX, EY, EZ - 4395007.97 4395007.97 1181426.56 
GY2,GXZ,GXY « 482009.43 482009.43 423769.42 
MUYZ,MUXZ,MUXY * .1274 .1274 .1370 


NUYZ,X / 

NUYZ,Y ; 

NUYZ,Z « 

. 0000 

.0000 

.0000 

NUXZ,X ; 

NUXZ,Y ; 

NUXZ,Z = 

. 0000 

.0000 

.0000 

NUXY,X ; 

NUXY, Y ; 

NUXY,Z = 

. 0000 

.0000 

.0000 
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